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Abstract 



The decays of massive gravitinos into neutralino dark matter particles and Standard Model sec- 
ondaries during or after Big-Bang nucleosynthesis (BBN) may alter the primordial light-element 
abundances. We present here details of a new suite of codes for evaluating such effects, including a 
new treatment based on PYTHIA of the evolution of showers induced by hadronic decays of massive, 
unstable particles such as a gravitino. We present several sets of results obtained using these codes, 
including general constraints on the possible lifetime and abundance of an unstable particle decay- 
ing into neutralino dark matter under various hypotheses for its decay mechanism. We also develop 
an analytical treatment of non-thermal hadron propagation in the early universe, and use this to 
derive analytical estimates for light-element production and in turn on decaying particle lifetimes 
and abundances, which confirm our numerical results and illuminate the underlying physics. We 
then consider specifically the case of an unstable massive gravitino within the constrained minimal 
supersymmetric extension of the Standard Model (CMSSM). We present upper limits on its possible 
primordial abundance before decay for different possible gravitino masses, with CMSSM parameters 
along strips where the lightest neutralino provides all the astrophysical cold dark matter density. 
We do not find any CMSSM solution to the cosmological 7 Li problem for small m 3 / 2 . Discounting 
this, for m-i/2 ~ 500 GeV and tan (3 = 10 the other light-element abundances impose an upper limit 
m 3 / 2 n 3 /2/ny < 3 x 10~ 12 GeV to < 2 x 10~ 13 GeV for m 3/2 = 250 GeV to 1 TeV, which is similar 
in both the coannihilation and focus-point strips and somewhat weaker for tan (3 = 50, particularly 
for larger m l / 2 - The constraints also weaken in general for larger m 3 / 2 , and for m 3 / 2 > 3 TeV we 
find a narrow range of m 3 / 2 n 3 / 2 /n 7 , at values which increase with m 3/ / 2 , where the 7 Li abundance 
is marginally compatible with the other light-element abundances. 
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1 Introduction 



Many extensions of the Standard Model of particle physics predict the existence of massive, 
unstable or metastable particles. If these have lifetimes longer than about a second, their 
decay products may affect the abundances of light nuclei, such as D, 4 He and 7 Li, produced 
during Big-Bang nucleosynthesis (BBN). If the particles have lifetimes exceeding about 100 
seconds, they decay after BBN, but their decay products may re-process the abundances 
produced during BBN. Both of these possibilities are strongly constrained by the degree 
of concordance between abundance predictions of D and 4 He and the observational deter- 
mination of these abundances, if the baryon density is within the range determined from 
measurements of fluctuations in the cosmic microwave background radiation (CMB) [1]. 
This concordance is not only an important vindication of BBN calculations [2-6] , but also 
sets severe limits on the possible abundance of any massive, metastable particle as a func- 
tion of its lifetime. These limits have been the subject of many previous studies that have 
modeled both the electromagnetic and hadronic components of the showers induced by the 
decays of such heavy particles [7]- [28]. 

We revisit here these constraints using a new suite of codes to model the decay spectra, 
the evolution of the showers that they induce in the electromagnetic plasma that filled the 
universe during the epoch between BBN and the decoupling of the CMB, and their potential 
for either destroying or creating light nuclear species [29]. Studying the decays of heavy 
particles during BBN requires a code that calculates the direct impact of electromagnetic 
and/or hadronic showers during nucleosynthesis, rather than just the effects on abundances 
that had previously been fixed independently by BBN. Our earlier code [20] considered only 
electromagnetic decays and did not provide constraints for particle lifetimes < 10 4 s. Here we 
describe a new suite of codes that includes hadronic decays and is applicable also to particles 
with lifetimes < 10 4 s. 

The effects included in this suite of codes are potentially relevant to any extension of the 
Standard Model that predicts the existence of massive unstable or metastable particles. One 
example is supersymmetry (SUSY). Many supersymmetric models incorporate a discrete Zi 
symmetry known as R parity. If R parity is not maintained, the lightest supersymmetric 
particle (LSP) is unstable, and the consistency of BBN with astrophysical observations 
and the density of baryons inferred from the CMB imposes important constraints on the 
mechanism and amount of R violation [30]. On the other hand, if R parity is conserved, 
the LSP is stable, but there is in general a heavier metastable supersymmetric particle. The 
possibility studied most frequently has been that the LSP is the lightest neutralino x [31]- 
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In this case, the gravitino would be metastable, and the BBN constraints on its decays 
impose important limits on the primordial production of gravitinos that are potentially 
problematic for some inflationary cosmological scenarios [32] - [40]. However, the neutralino 
is not the only candidate for the LSP, and another generic possibility is the gravitino itself. 
In this case, the gravitino itself would be stable, but the next-to-lightest supersymmetric 
particle (NLSP) would in general be metastable, and its decays would be subject to the 
BBN constraints [41-43]. 

Our goals in this paper are threefold. The first is to document our new interlocking suite 
of codes for decay showers in the early universe and their impact on BBN; our codes include 
the effects of both hadronic and electromagnetic showers, and are applicable to a large range 
of unstable particle lifetimes. The second goal is to set out the constraints obtained using 
this code on the lifetime and abundance of a generic unstable or metastable particle, under 
various plausible hypotheses for its dominant decay modes. The third goal is to explore 
in more detail the consequences of the shower effects for supersymmetric models with a 
neutralino LSP. In the context of such models, we demonstrate how our new BBN codes 
may be used to obtain constraints on the gravitino abundance as a function of its mass. 

In previous analyses of the constrained minimal supersymmetric extension of the Stan- 
dard Model (CMSSM) with a neutralino LSP, cosmological constraints on the decays of the 
metastable gravitino have usually been ignored. Essentially, it has been assumed implicitly 
that the gravitino is so heavy, and hence its lifetime is so short, that its decays take place 
so early in the history of the universe that they do not affect light-element abundances. 
However, from a theoretical point of view, a large hierarchy between the gravitino and LSP 
masses may appear unlikely. Indeed, in minimal supergravity (mSUGRA) and related mod- 
els [44], it is usually the case that the gravitino is not much heavier than the LSP, if it is not 
the LSP itself. Therefore, it is important to analyze the case when the decays of gravitinos do 
affect the light-element abundances, and evaluate the resulting constraint on the primordial 
gravitino abundance as a function of its mass. 

We exemplify these constraints in the specific case of the CMSSM, in which supersymmetry- 
breaking gaugino and scalar masses mi/ 2 and m are each assumed to be universal [45,46]. In 
this case, when the LSP is the lightest neutralino, its relic density 0.0975 < f2 C DM^ 2 < 0.1223 
as inferred from astrophysical and cosmological measurements [1] constrains (mi/2, m ) to lie 
along narrow "WMAP strips" for given fixed values of the other CMSSM parameters [45,47]. 
Moreover, the gravitino lifetime and decay modes are known as functions of its mass and the 
CMSSM parameters. Thus, we are able to set firm upper limits on the possible abundance of 
the gravitino as a function of its mass and mi/ 2 . This upper limit on the gravitino abundance 
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may, in some circumstances, constrain severely the maximum temperature reached in the 
universe, e.g., following an early inflationary epoch [32] - [40]. 

The layout of this paper is as follows. In Section [2] we summarize the cosmological data 
that we use in our analysis. Then, in Section [3] we explain our treatment of the evolution 
of the hadronic and electromagnetic showers initiated by the decay of a generic metastable 
particle into the cosmic plasma, and their effects on the primordial light-element abundances. 
Appendix [A] summarizes our treatment of hadronic spectra, and describes the propagation 
of non-thermalized particles in the cosmic plasma. In Section H] we turn to the specific case 
in which the gravitino is the metastable particle, relegating some technical details to Ap- 
pendix [HI In Section [5] we discuss the resulting constraints on the abundance of a metastable 
neutral particle, under plausible assumptions on its dominant decay modes. We consider con- 
straints as a function of a generic particle lifetime, as well as constraints on the gravitino 
abundance in specific CMSSM scenarios with a neutralino LSP. For m 3 / 2 < 1 TeV, we find 
no parameter choices where the 7 Li abundance can be reconciled with the other light-element 
abundances. However, for m 3 / 2 > 3 TeV we find narrow ranges of the gravitino abundance 
that may reconcile marginally the 7 Li and other light-element abundances. Finally, in Section 
|6] we summarize our conclusions and present some prospects for future work. 

2 Cosmological Data 

2.1 Light-Element Abundances 

The abundances of the light elements D and 4 He predicted by BBN theory agree quite well 
with the values determined by observation, if the baryon-to-photon ratio rj is that inferred 
from CMB measurements [1]. This concordance provides the basis for the constraints on 
metastable particles to be discussed in this paper. However, there is known to be an issue 
regarding the abundance of 7 Li, which we discuss below. 

Deuterium provides a powerful constraint, as it is very sensitive to the baryon content 
in the universe, and thus offers by itself a measure of rj. Local deuterium that is measured 
in the solar neighborhood in the interstellar medium provides a strong lower limit on the 
cosmological abundance, given that astrophysical effects destroy more deuterium than they 
produce [48]. Recent observations by FUSE show a wide dispersion in the deuterium abun- 
dance in local gas seen via its absorption, (D/H) local gas = (0.5 — 2.2) x 10~ 5 [49]. This surpris- 
ingly large spread, taken together with the positive correlation of D /H with temperature and 
metallicity along various sightlines, led [49] to suggest that deuterium may suffer significant 
and preferential depletion onto dust grains. In this case the true local interstellar D/H value 
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would lie at the upper limit of the observed values, giving (D/H) ISM > (2.31 ± 0.24) x 1CT 5 . 
However, extracting a primordial deuterium value requires a Galactic chemical evolution 
model (e.g., [50]), whose model dependences yield uncertainties in the determination of the 
primordial deuterium abundance. 

A high-redshift, metal-poor system is free of many such model uncertainties. In partic- 
ular, observing the absorption features of quasar light due to a dense high-redshift cloud 
can be used to determine the primordial deuterium abundance. Several studies yield re- 
sults [51-57] that are broadly consistent, but with a reduced xl — 2-95, suggesting that the 
uncertainties have been underestimated. Taking a weighted mean, we find a "world average" 
range: 



where the error bar has been inflated by a factor of \fxl ■ Because it is likely that systematic 

errors dominate, an even more conservative approach would be to use the sample variance 
of the D/H data [2]; this would give a(D/H) = (±0.53) x 10~ 5 , i.e., a significantly higher 
error. In this study we adopt 3.2 x 10~ 5 as our fiducial upper limit on the D/H abundance, 
but we also illustrate the effect of significant variations in D/H around this value. 

Since 3 He is also quite sensitive to the baryon density, one might hope that it too could 
be used as a baryometer. Observations of HII regions in our own Galaxy yield values of the 
3 He/H ratio that are compatible with calculations of the primordial value [58,59]. However, 
the extrapolation from the observations to a primordial abundance is complicated by the 
unknown chemical evolution of 3 He. Indeed, one does not even know whether 3 He/H is 
increasing or decreasing with cosmic time. Thus, a primordial extrapolation yields only an 
order-of-magnitude range of allowable values of 3 He/H [60]. However, if we use the results 
of [48, 61] that the deuterium abundance is always decreasing with time, and assume that 
3 He changes relatively slowly, we can adopt their ratio: 



as a conservative constraint on the primordial 3 He/D ratio. 

HII regions also yield observations of 4 He. However, measurements of extra-galactic and 
metal-poor sources yield a priori more reliable estimates of the primordial 4 He abundance. 
Various evaluations of the same data yield differing results ranging from Y p = 0.234 [62] to 
Y p = 0.244 [63]. These differences point to the dominance of systematic uncertainties [64] 
over statistical errors. The allowed ranges of systematic shifts are discussed in [65], but their 
analysis does not provide a full estimate for the systematic error. Subsequently, another 
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analysis of the data was performed [66] , making more realistic assumptions than those used 
in previous analyses. This analysis found a systematically higher Y p , with significantly 
increased errors, suggesting that previous analyses had underestimated their systematics: 



Y p = 0.249 ± 0.009. 



(3) 



The inclusion of one of these systematics has been explored in [67], which finds results similar 
to [66]. In our subsequent analysis, we adopt the lower limit Y p > 0.240. 

Lithium is seen in the atmospheres of the most primitive, metal-poor stars in the stellar 
halo of our Galaxy (extreme Population II). Some > 100 such stars show a plateau [68] 
in (elemental) lithium versus metallicity, with a small scatter consistent with observational 
uncertainties. This insensitivity of Li/H to (proto)-Galactic metal content and thus stellar 
nucleosynthesis indicates that lithium in these stars has a primordial or at least pre-Galactic 
origin. An analysis [69] of field halo stars gives a plateau abundance of 



where the errors represent a 95% CL and include both statistical and systematic uncer- 
tainties. On the other hand, studies of globular cluster stars generally find higher lithium 
abundances. The measurement [70] of 



is consistent with results for other systems: Li/H = (1.91 ± 0.44) x 10~ 10 [71]; Li/H = 
(1.69 ±0.27) x 10- 10 [72]; and Li/H = (2.29 ±0.94) x 10~ 10 [73]. Both sets of abundances are 
substantially lower than the standard BBN prediction of 7 Li based on the WMAP estimate 
of i] [1], by factors of ~ 2.4 — 4.3, or 4 — 5a [6]. This discrepancy is known as the "lithium 
problem" or more specifically the " 7 Li problem." In our later analysis, we adopt 7 Li/H < 
2.75 x 10~ 10 as our fiducial upper limit on the cosmological 7 Li abundance. 

The existence of the 7 Li problem, and the nature of its solution, both have a direct 
impact on our analysis. It is possible that the problem points to new physics, in particular 
if the observations and standard theory are both correct with accurate error estimates. If 
SUSY were to lead to a solution of the 7 Li problem, this would tie together a wide array 
of particle astrophysics experiments and observations. Indeed, it has been suggested [21,22] 
that the 7 Li problem could be solved by hadronic decays of a metastable neutral particle X 
with lifetime t x ~ 10 3 sec. We will directly address and update this issue below. On the 
other hand, it remains possible that the standard BBN light-element abundance predictions 




(4) 




(5) 
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remain correct, i.e., unperturbed by any dark matter interactions and/or decays. Rather, 
it could be that the 7 Li problem instead reflects astronomical observational systematics in 
recovering a Li/H abundance from stellar spectra, though a recent study [74] of the effective 
temperature of metal poor stars confirmed the use of relatively low temperatures and a Li 
abundance in the range Li/H = (1.3 — 1.4 ± 0.2) x 10~ 10 . The 7 Li problem may also reflect 
astrophysical systematics due to Li depletion via circulation and nuclear burning over the 
> 10 Gyr lifetime of metal-poor stars [75], though the lack of star-to-star scatter in Li/H 
suggests the depletion could be negligibly small [69]. If there is depletion, for our problem 
the correct procedure would be to use the initial, undepleted Li/H abundances. However, 
since these cannot be determined accurately, we do not adopt a quantitative constraint for 
this possibility. Instead, we consider the limits on unstable particles when the 7 Li is assumed 
to be solved and thus the 7 Li constraints are ignored. 

Recent work has claimed to resolve lithium isotope information in halo stars [76,77], and 
indicates the presence of 6 Li. In fact, the isotopic ratio has a roughly constant value of 

6 Li 

y— «0.05. (6) 

This has several immediate implications. First, we see that most of the lithium is in the from 
of 7 Li, so that the total elemental lithium abundance (as in eqs. (j3J) and is mainly 7 Li: 
Li/H w 7 Li/H. Secondly, given that Li/H in these stars has a plateau, the constancy of the 
isotopic ratio implies that both isotopes have a plateau. However, the inferred 6 Li plateau 
abundance is about 1000 times the 6 Li abundance predicted by standard BBN [78,79]. 

Extracting isotopic information from thermally broadened absorption lines is extremely 
challenging, and the existence of a 6 Li plateau has been questioned [80]. If real, the 6 Li 
plateau abundance cannot be explained by conventional Galactic cosmic-ray nucleosynthesis 
[79,81,82]. Therefore, one is led to think that there is a " 6 Li problem" in addition to that 
for 7 Li. 

It is possible to arrange decaying-particle scenarios in which 6 Li is produced at the plateau 
level with some destruction of 7 Li, thus fixing both lithium problems simultaneously [26,29, 
83-86], a point to which we will return in detail below. However, even granting the existence 
of the 6 Li plateau, the 6 Li problem really only requires a pre-Galactic source, which need 
not arise prior to the first star formation. For example, the 6 Li problem might be explained 
by cosmological cosmic-ray nucleosynthesis due to cosmic rays produced at the epoch of 
structure formation [87]. For our purposes, therefore, we simply interpret eq. (jSJ) as a firm 
upper limit on any primordial 6 Li [77,80]. 
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2.2 Baryon Density 



The baryon-to-photon ratio rj = Ub/ti^ = 10~ 10 ?7io is related to the present baryon density 
parameter Qb = Pb/ Pcrit by ^b^ 2 = ?7io/274. Adopting the value of VLbIi 2 indicated by the 
WMAP 5-year CMB data [1] gives the following estimate of the baryon-to-photon ratio: 



This is the default value we assume later in our BBN analysis. The central value and error 
range here correspond to WMAP data combined with large-scale structure information, in a 
framework for which the primordial power spectrum is a simple power law with fixed index. 
Similar but slightly different values would result from, e.g., a running spectral index. 

3 Hadronic Decays during Primordial Nucleosynthesis 

Particle decays during BBN generally have two main effects [7,31]. First, they change the 
cosmic expansion rate due to the injection of additional relativistic species. This effect is 
model-independent, in that it is insensitive to the details of the decays, beyond the as- 
sumption that the daughter particles are relativistic [10]. Secondly, they introduce new, 
non-thermal decay products-possibly electromagnetic and/or hadronic-which can interact 
with the background thermal nuclei and change the final light-element abundances. The 
branching ratios and spectra of hadronic and electromagnetic decay particles and energies 
are model-dependent, and are further discussed below in Section HI While both effects occur 
in principle, the expansion effects are negligible for the decay parameters of practical interest 
to us, and the non-thermal interaction effects are the dominant perturbations to the abun- 
dances. Consequently, we focus on these effects, and henceforth neglect the perturbation 
to the expansion rate. If the decaying particles are electrically charged and thus can form 
bound states with nuclei, additional effects arise [29,83-86,88-96]. This is not the case in 
the scenarios considered in the present paper, but we will revisit this subject in future work. 

To fix notation: we consider decays of some generic heavy particle X, which we will 
eventually specialize to the case of the gravitino. The particle has lifetime Tx and decay rate 
Fx = t^ 1 - We can quantify the X abundance either as the number of decaying particles per 
background baryon, 



mo 



6.23 ±0.17. 



(7) 
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n x v n B 1 
= — = Y x — « -rjYx, 
s s 7 



(9) 
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where s ~ 7n 7 is the entropy density, and the factor of ~ 7 is appropriate after e annihila- 
tion. We note that, in the absence of X decays, conservation of baryon number guarantees 
the constancy of Yx, while yx changes during BBN due to the usual photon heating from 

annihilation. Due to decays, we have Yx(t) = Y^ lt e~ t l Tx \ our constraints will be on 
the initial, pre-decay abundance Y^ mt , and this should be understood hereafter whenever we 
write Yx- The only exception in our notation is that we write Y p for the 4 He mass fraction 
and this should not be confused with the ratio of number densities as defined in eq. (jHJ). 

We are interested in both hadronic and electromagnetic decays of the heavy particles X. 
We denote the electromagnetic branching ratio of X by .Bem = ^ x^~em/^ x- The abundance 
perturbations from electromagnetic decays, and thus the constraints on such decay modes, 
scale with the product of 5em and the decay energy release per photon: 

. m x n x v , v 

(x = = m x Y x V: (10) 

where the X abundance is evaluated prior to decay. As emphasized particularly by [22], 
electromagnetic decays inevitably accompany hadronic decays, and so both sets of decay 
products and interactions need to be included. Electromagnetic decays were discussed in 
detail in [29]; we include those processes here, incorporating the treatment described in [29]. 
However, the analysis there assumed that the electromagnetic decays occurred entirely after 
the usual BBN thermonuclear reactions have run to completion, i.e., the decays were treated 
as a "post-processing" modification after the usual light-element abundances had been es- 
tablished. Here, we supplement the previous treatment by including electromagnetic decay 
effects consistently throughout BBN. 

We quantify generic decays X — > h of a particle X into a hadronic species h as follows. We 
write the h production rate per unit volume, and daughter kinetic energy e as q7i(e). The total 
X decay rate per baryon is YxYx and the total X decay rate per volume is qx,tat — ^xYxTIb- 
It will be convenient to isolate the h production per X decay as Qh{e) = qh(e)/qx, tot, so that 
we have 

q h (e) = T x Y x n B Q h (e). (11) 

We refer to Qh{t) as the spectrum of X decays, which gives the number of h particles 
produced per X decay per energy interval. The integral Bh = / Qh( e ) de gives the total 
number of h per X decay, and thus represents a multiplicity-weighted hadronic branching 
ratio. Gravitino decays are described in detail in the next Section. 
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3.1 Interactions with Background Nuclides 

The dominant effect of hadronic decays on BBN is the addition of new interactions between 
hadronic shower particles and background nuclides. These act to alter the evolution and final 
values of the light-element abundances, as follows. For each light nuclide £, the abundance 
per background baryon, Y e = ni/n B , changes according to 

d t Y E = (d t Y e ) SBBN + (d t Y e ) EU + («9^) had , (12) 

where (<9^)sbbn denotes the usual rate of change of £ in standard BBN due to thermonu- 
clear reactions. The second term on the right-hand side of (Tl2l accounts for non-thermal 
electromagnetic interactions due to X decays, either directly from the decays of X to pho- 
tons or leptons, or through electromagnetic secondaries in the hadronic showers. These are 
treated as in [20], but are not dominant when hadronic branchings are significant. The final 
term on the right-hand side of ( fl2|) represents the non-thermal hadronic interactions, which 
are a major focus of this paper. 

The inclusion of hadronic decays in BBN thus requires a detailed evaluation of the rates 
for such interactions. For each background light nuclide species £, we can write the hadronic 
contributions to d t Ye as 

(<9^)had = — r^ineiY^ + ^^T hb ^gY b . (13) 

hb 

The first term accounts for £ destruction by hadro-dissociation, where r^i ne i is the total rate 
(per unit of the species £) of all inelastic interactions of shower particles with £. The second 
term accounts for production due to hadro-dissociation of heavier background species (e.g., 
deuteron production via helium erosion p s howe r abg — ► d + ■ • •). The sum of inelastic rates 
^hb-*t producing £ runs over shower species h and background targets b. In the particular 
case of lithium isotopes, production occurs via the interaction of energetic (i.e., non-thermal) 
mass-3 dissociation products with background 4 He, e.g., 3 He + a — > 6,7 Li + • • • . 

The non-thermal reaction rates on light elements are themselves set by the abundance 
and energy distribution of non-thermal particles, which we quantify as follows. Consider 
a hadronic (non-thermal projectile) species h, with an energy spectrum ^|-(e, t) (particle 
number per unit volume per unit energy interval) and total number density rih = J ^-de. 
It is convenient to define the energy spectrum of non-thermal h as = which is 

normalized such that f N h (e) de = n h /n x = n h /Y x n-Q measures the number of propagated 
non-thermal h particles per X . 
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The rate for non-thermal production of species t due to hb — > £ is 

T hb ^t = Y x n B J N h (e,t) a hb ^ e (e) f rc i(e) de « Y x n B J N h (e,t) [a hb ^ e v] (e) de, (14) 

where ^^(e) is the cross section for this process. Here f rc i(e) ~ Vh(e) is the relative 
velocity between the projectile and the background target, which in practice amounts to the 
projectile velocity. In a similar way, the inelastic-loss rate is just a sum over all inelastic 
channels removing £, i.e., where the final state / does not include £: 

r^mci ~ Y x n B ^ / N h (e,t) [a M ^f v ] (e) de. (15) 
/ J 

The rates Yi thus depend on the shower development and evolution of Nh(e,t) in the back- 
ground environment. A major effort of this paper is to calculate these spectra and their 
evolution. 



3.2 Hadronic Showers 

We wish to follow the evolution of the non-thermal spectra Nh over the multiple genera- 
tions produced by the hadronic X decays. In the context of BBN, this problem of shower 
development has been approached via Monte Carlo computation of the multiple generations 
of shower particles [18,19,21,22]. The final particle spectrum is obtained via iterating an 
initial decay spectrum, accounting for both energy losses as well as the energy distribution 
of collision products. 

Here we introduce an alternative, but equivalent, approach based on a "cascade equation" 
treatment, that follows the well-studied treatment of the development of hadronic showers 
due to cosmic-ray interactions in the atmosphere. As we show, this approach offers new 
physical insight as well as computational advantages. 

Once a hadron h is produced by X decay, it loses energy, interacts with background nuclei, 
and possibly decays if it is unstable (e.g., n, 7r ± ) . The spectrum of h evolves according to 
the energy-space "propagation equation" : 

d t N h (e) = J fc (e) - T h (e)N h (e) - d e [b h (e)N h (e)] , (16) 

where Jh is the sum of all source terms, is the sum of all sink terms, and bh is the 
energy-loss rate of particle-conserving processes. Note that inelastic interactions lead to 
secondary production of h, sometimes by other non-thermal species h! . For this reason, 
eq. ([TBI ultimately represents a set of coupled, integro-differential equations. 
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Details of eq. (Fl6|) and its numerical solution appear in Appendix |A] In the next section, 
we summarize the key physics via an analytical description informed by the full numerical 
results. 

3.3 Non-Thermalized Particle Spectra and Interactions: Analyti- 
cal Model 

To develop a physical intuition for the full numerical results used later, we first present 
a simplified analytical treatment of the propagation and interactions of non-thermalized 
hadrons. 

A crucial feature of eq. (1T61) is that it contains both source and sink terms for h. The rates 
for energy loss in the cosmic plasma, and for scattering on the background nuclei, are always 
faster than the cosmic expansion rate H. Thus, redshift and expansion effects may safely be 
ignored in non-thermal particle propagation, and have been neglected in eq. ( 1T61) . Moreover, 
the energy-loss and scattering rates are generally faster than the non-thermal source rate 

— 1/tx set by the X decay time. Thus eq. ffl6l) is self-regulating, with tracking a 
quasi-static equilibrium set by the balance of source and sink terms. That is, at any cosmic 
epoch we have dfNh(e) ~ 0. This provides a crucial simplification and reduces the problem 
to a series of ordinary integro-differential equations, which one may solve at each epoch. 

We thus begin with a simplified version of the full propagation equation, which neglects 
elastic and inelastic source terms: 

= -d e (b h N h ) - (I\ sc + T p )N h + T x Q h , (17) 

where bh = —deh/dt is the total energy-loss rate, and I\ sc is the total rate for inelastic and 
elastic scattering. For neutrons, Tp = l/7„r n is the rate for the beta decay of a free neutron 
with Lorentz factor j n ; for protons Tp = 0. 

Equation (TlTl) is formally identical to the "leaky box" equation for the propagation of 
cosmic rays in our Galaxy. Indeed, our problem can be regarded as the injection of cosmic 
rays in the early universe plasma. The solutions to eq. ( TlTl) have been well-studied for the 
case of Galactic cosmic rays [97], and we adapt the solutions for our problem. We focus on 
the propagation of the primary species, i.e., the p and n, which have sources in the X decays 
themselves; in this case the decay term dominates the sources. Our full numerical treatment 
also includes primary tt^, but these are important only for short Tx ^ 1 sec. 

On physical grounds, we expect that the hadronic spectra at any given time should 
scale linearly with the X decay rate Tx- One can infer this formally from eq. (Tl6|) . as outlined 
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in Appendix [A] Thus we have the scaling Nh oc Fx, and so the ratio Nh/Tx is independent 
of the decay rate Tx- 

If we ignore secondary production of h, the exact solution to eq. ( fl7l) is known and 
appears in Appendix |A] However, for our purposes it is useful to focus on two limiting cases. 
Generally, one of the sink terms is dominant and controls the resulting equilibrium. If the 
scattering term dominates the energy- loss term, FN ^> d e (bN) we refer to this as the thin- 
target limit. Physically, h particles interact catastrophically before they lose energy, i.e., the 
cosmic plasma is "thin" against stopping. The reverse case, TN <C d e (bN), constitutes the 
thick-target limit, where particles lose energy before interacting. We can express the ratio of 
the two terms as d e (bN)/TN = (b/eT)dln(bN)/dlne. We see that, as long as the logarithmic 
term is slowly varying, the dominant loss term is set by a comparison of scattering timescale 
r _1 with energy-loss timescale e/b: 

T sc > - thin target, (18) 

r sc < - =^ thick target. (19) 

In the thin-target limit, scattering losses dominate, so that eq. (TTTj) reduces to —T h)SC N h + 
TxQh — and the solution becomes algebraic 

n,M - 1* ^ = r. Y Y Qh(e) .. . (20) 

Here the reaction on background species b dominates the scattering losses. Note the inverse 
scaling with baryon density n&. Note also that we have implicitly assumed rg <S T sc . In 
the case of neutrons, once Tp > r inol (which is indeed well into the thin-target regime), 
then neutrons decay before they interact. This occurs when r ine i n ~ n B v n a ine i < Tp, which 
corresponds to T < 0.4 keV 7^ 1/3 (w„/0.1c)- 1 / 3 and t > 8 x 10 6 s ^ /3 {v n /0.1c) 2 / 3 . At longer 
times, non-thermal protons and, more importantly, electromagnetic cascades, dominate light- 
element production. 

Turning now to the thick-target limit, when energy losses dominate the sinks we have 
~ -d e (bN h ) + T x Q h , and thus 

N h (e)^r x ®^, (21) 
6(e) 

where Qhi> e) = / Qh( e ') dd is the integral source spectrum above e. Taking b = b con i for 
Coulomb losses, we have b ~ 47ra 2 h 2 c 2 Y e ±n-B/m e v , or 

de/dR = b/p B v ~ 47ra 2 h 2 c 2 Y e ±/m e m p v 2 ~ 1.5MeV/(gcnT 2 ) F e ±(100MeV/e). (22) 
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Crucially, the electron/positron abundance Y e ± per baryon is enormous before pair annihi- 
lation, so that Y e ± ~ l/rj ~ 10 9 at T > 0.2 MeV, and even afterwards the pairs dominate 
the baryons until T ~ m e /25 ~ 0.02 MeV. Thus we write 

m p nBV{e) deh/dR{e) 

where we use the energy loss per grammage de/dR = b/psv, which is independent of back- 
ground density. Here again we see that the spectrum scales inversely with background den- 
sity, but also has an inverse dependence on the strongly- varying electron/positron abundance 
Y e ±, contained in de/dR. 

3.4 Non-thermal Reaction Rates 

Consider a reaction hb — > £ on a light background species b which produces a daughter 
species £, e.g., na — > d H . The rate for this reaction is 

T hb ^e = J <f>h(e) Vhb^e(t) de (24) 

per background target b. Here the (angle-integrated) flux of non-thermal h e n, p is ^^(e) = 
nxNh(e)v(e) = Yxn-QNh{e)v{e). Recalling that iVj, oc 1/ne, we see that the non-thermal flux 
4>h oc nsNh and thus the reaction rate Thb^e are both independent of the background baryon 
density. Also, since Nh oc Tx, we have Thb^e oc TxYx- 

We now examine the reaction rate for the two limiting cases of propagation. Using Nh 
in the two limits, we have 



hb^i — 



f Qh(e) de thin target, 

r x lx r £th Q,(> e) a -^f^ de thick target, 



where the integral begins at the threshold energy e t h (if any). Here we take the interaction 
hk — > inelastic of the non-thermal particle ft, with background nuclide to dominate the 
inelastic losses. 

With these reaction rates in hand, we illustrate their importance by estimating below 
the perturbation on light-element abundances that arises in the CMSSM with a massive 
gravitino. Before doing this, we must first have an understanding of the gravitino decays 
and the resulting hadronic spectra. 
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4 Gravitino Decays 



We work in the context of the CMSSM [46] which is described by four parameters: universal 
gaugino, scalar, and trilinear masses, mi/2,mo, Ao, and the ratio of the two Higgs vacuum 
expectation values, tan (3, along with the sign of the Higgs mixing parameter, [i. Motivated 
by g^, — 2 and h — > 57, we restrict our attention to > and, for simplicity, we consider 
only A = 0. We consider scenarios in which the lightest neutralino is the LSP, but the 
gravitino is not assumed necessarily to be the NLSP. For fixed values of tan/3, the regions 
of parameter space for which the relic density of neutralinos is computed to lie within the 
range 0.0975 < f^DM^ 2 < 0.1223 determined by WMAP [1] and other observations for 
cold dark matter form narrow strips that foliate the (mi/2,^0) plane [45,47]. These strips 
correspond to regions where there are enhanced annihilation cross sections that reduce the 
neutralino relic density to acceptable values. These strips occur when the neutralino is nearly 
degenerate with some other supersymmetric particle such as the partner of the tau lepton 
(coannihilation strip), when the neutrino is close to half the mass of the heavy Higgs scalar 
and pseudoscalar boson so that rapid s-channel annihilation occurs (the funnel region), or 
at large values of tuq when the renormalization-group evolution drives the value of /1 to low 
values so that the neutralino acquires a significant Higgsino component and new final-state 
channels become important (the focus-point strip). Our analysis is based on the WMAP 
strips for two representative values of tan (3 = 10, 50. For tan (3 = 10 one strip follows the 
coannihilation corridor, and for tan (3 = 50 this strip also includes the funnel at larger values 
of mi/2- We also consider the focus-point strips for both values of tan/3. 

In the scenario under study the gravitino decays into lighter sparticles, including the 
neutralino LSP but also others in general. We take into account all the dominant decay 
channels of the gravitino into (s)particles, including the complete set of two-body decays of 
the gravitino into Standard Model particles and their spartners. These fall into the following 
main categories: G -f / /, G -> x + W~(H~), G -> $7(Z), G -> x° Hf and G -> gg. 
Analytical expressions for these amplitudes can be found in Appendix [B] In addition, we 
include the dominant three-body decays G — > xl 7* — ► xl QQi G — > xl W + W~ , which are also 
discussed in Appendix El In principle, one should also include qq pair production through the 
virtual Z-boson channel G — > Xi Z* — > xl, QQ [24] and the corresponding interference term. 
However, this process is suppressed by a factor of M| with respect to G —>■ xl 1* ~ > Xi QQi an d 
the interference term is also suppressed by Mf . These contributions are therefore not very 
important, and we drop these amplitudes in our calculation. We also drop the corresponding 
amplitudes for Higgs and squark exchange. We calculate the lifetime of the gravitino by first 
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calculating the partial widths of its dominant relevant decay channels, and then summing 
them. Typical results are shown in Fig. [T], where we see that the lifetimes for tan f3 = 50 
(right panels) are typically longer than those for tan/3 = 10 (left panels), particularly for 
larger m^j. On the other hand, the lifetimes in the coannihilation/funnel regions (top 
panels) are quite similar to those in the focus-point regions (bottom panels). 

In order to calculate the resulting electromagnetic (EM) and hadronic (HD) spectra, we 
first calculate the EM and HD decay spectra of the different gravitino decay modes, and 
then weight them by the corresponding branching ratios. The decay products that yield EM 
energy obviously include directly-produced photons, but also indirectly-produced photons, 
charged leptons (electrons and muons) and neutral pions (vr°), which are produced via the 
secondary decays of unstable heavy particles such as gauge and Higgs bosons. Hadrons 
(nucleons and mesons such as the K*[, K and 7r ± ) are also usually produced through the 
secondary decays of heavy particles, as well as (for the mesons) via the decays of the heavy r 
lepton. It is important to note that strange mesons and neutral pions decay before interacting 
with the hadronic background [19,41,42]. Hence they are relevant to BBN processes and to 
our analysis only via their decays into photons and charged leptons, which contribute to the 
EM component of the decay showers. Therefore, the HD injections that concern us are those 
that produce nucleons, via the decays of heavy particles such as gauge and Higgs bosons and 
quark- ant iquark pairs and to a lesser extent charged pions. 

After calculating the partial decay widths and branching ratios, we then employ the 
PYTHIA event generator [98] to model both the EM and the HD decays of the direct products 
of the gravitino decays. We first generate a sufficient number of spectra for the secondary 
decays of the gauge and Higgs bosons and the quark pairs. Then, we perform fits to obtain 
the relation between the injected energy of the decaying particle, that we call Ei n j, and the 
quantity that characterizes the hadronic spectrum, namely Qh, the number of produced nu- 
cleons as a function of the nucleon energy, for various values of the E in j. We have performed 
fits that cover the range 200 GeV < Ei n j < 20 TeV. In Table [1] we present a representative 
set of fitting parameters for E in j = 1000 GeV. Then, using the equation 

Q h (e h ) = ^-Q h (x), (26) 

-G/ in j X 

the spectrum distribution Qh{^h) is computed using the function Qh(x) = Qh(€h)deh/dx, 
with x = a/c| — ml/ E in j. Some typical spectra for gravitino decays into protons are shown 
in Fig. [2J These spectra and the fraction of the energy of the decaying particle that is 
injected as EM energy are then used to calculate the light-element abundances. We stress 
that this procedure is repeated separately for each point sampled in the supersymmetric 
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Figure 1: The gravitino lifetime for representative values of 1713/2 as a function of m 1( / 2 
for tan (3 = 10 (left) and tan/? = 50 (right) along WMAP strips, in the coannihilation and 
funnel region (top) and the focus-point region (below). 



parameter space. That is, given a set of parameters m , m^, A , tan/3, sgn(/i), and m 3 / 2 , 
after determining the sparticle spectrum, all of the relevant branching fractions are computed, 
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and the hadronic spectra and the injected EM energy determined case by case. Thus, in 
our analysis, Q p (Q n ) and hence the total number of protons (neutrons) per gravitino decay, 
Bp (B n ), varies between different points in the parameter space, and Fig. [3] illustrates this 
variation in B p (and B n ) across the supersymmetric parameter space we sample. 

Table 1: The values of the parameters A, B, C , D and E used in fitting PYTHIA hadronic 
decay spectra. We use the function Y = log 10 (Q(a;)) parametrized via Y = AX 4 " + B X 3 + 
C X 2 + D X + E , where X = log 10 (x), and the scaling parameter x is defined as x = 
rn^/Einj. The values of this table correspond to E in j = 1 TeV. 



decaying particles 


inj. species 


A 


B 


c 


D 


E 


Z 


P 


0.000 


0.000 


-1.230 


-4.823 


-3.401 




n 


0.000 


0.000 


-1.281 


-4.917 


-3.418 


h 


p 


0.000 


0.000 


-1.011 


-4.302 


-3.071 




n 


0.000 


0.000 


-0.850 


-3.628 


-2.499 


H 


p 


0.000 


0.000 


-0.647 


-3.398 


-2.154 




n 


0.000 


0.000 


-0.644 


-3.414 


-2.230 


A 


p 


0.000 


0.000 


-0.633 


-3.356 


-2.133 




n 


0.000 


0.000 


-0.637 


-3.377 


-2.179 


qq 


p 


-0.065 


-0.551 


-2.007 


-4.326 


-2.602 




n 


-0.033 


-0.195 


-0.603 


-2.039 


-1.327 


w + w- 


p 


0.000 


0.000 


-1.146 


-5.180 


-4.036 




n 


0.000 


0.000 


-1.107 


-4.936 


-3.742 



The relative sizes of the gravitino partial widths and the locations of the various particle 
thresholds help us to understand the lifetime and the hadronic spectra curves. In general, 
the two-body channels G — > x° 7; f f] 9 9] X + W~ dominate the G decays. In particular, the 
decay to xl, 7 yields the bulk of the injected EM energy. However, the decays to sfermions, 
gluinos and charginos become of the same order of magnitude as the G — > channel 
whenever they are kinematically possible. Also, when the G is heavy enough to decay into 
a real Z boson, the channel G — > Xi Z is the dominant channel for producing HD injections. 
When kinematically allowed, G — > x + W~ and G — > g g are also important in producing 
HD injections. The Higgs boson channels are smaller by a few orders of magnitude (due 
to couplings and kinematics) and, in particular, in the large-mi/2 region, decays to heavy 
Higgs bosons (H, A) become kinematically accessible only for heavy G and are unimportant 
otherwise. 

Turning to the three-body channels, the decay through the virtual photon to a qq pair 
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Figure 2: Sample spectra for gravitino decays into protons. We plot values of the combi- 
nation e N p (e), which gives the particle number per logarithmic energy range, for different 
representative values of the gaugino mass parameter myz along the coannihilation strip for 
tan/? = 10, assuming the indicated values of the gravitino mass m 3 / 2 - 

can become comparable to the channel G — > Xi Z, injecting nucleons even in the kinematical 
region 777,3/2 < m x + Mz, where direct on-shell Z-boson production is not possible. Finally, 
we note that the partial width of the three-body decay G — > x° W + W~ is usually smaller by 
at least an order of magnitude relative to the dominant two-body decays, except when this 
three-body decay exhibits resonant behavior. For example, the subprocess G — > x + *W~ — > 
Xi W + W~ can lift the contribution of the %° W + W~ channel to the level of the dominant 
two-body decays in the threshold region where the chargino can be produced on-shell. 

With these observations in mind, one can understand the gravitino lifetime curves along 
the WMAP strips for tan (5 = 10 and 50 in Fig. [TJ We recall that the funnel region is only 
present in the CMSSM for large tan/3. Since the relation m x w 777,4/2 is realized at large 
mi/2, the WMAP strip extends to significantly higher values of mi/ 2 for tan/5 = 50 than 
for tan/5 = 10. In the latter case (upper left panel) , the WMAP strip shown consists only 
of a coannihilation region, which terminates around myz = 900 GeV. We notice that the 
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Figure 3: The number of nucleons per gravitino decay, as a function ofm\ii for tan (3 = 10 
(left) and tan/3 = 50 (right). The upper panels are for WMAP strips in the coannihilation 
and rapid- annihilation regions, and the lower panels are for WMAP strips in the focus-point 
regions. Solid curves: number B p of protons per decay. Broken curves: number B n of 
neutrons per decay. We see that generally B p « B„ to a good approximation. 
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gravitino lifetime is longer for tan (3 = 50 than for tan (3 = 10, for the same values of 7B3/2 and 
mi/2- Especially for the two upper panels, this is because for tan/3 = 50 the WMAP strip, 
in the coannihilation region and (particularly) in the Higgs rapid-annihilation funnel region, 
occurs at larger m and hence heavier squark and slepton masses, than in the tan/3 = 10 
case. This implies that the dominant two-body channels G —>■ f f are very suppressed or 
even closed for tan (3 = 50. Thus, the total gravitino decay width is smaller (and hence the 
lifetime longer) for tan j3 = 50 than for tan f3 = 10. 

The lower panels in Fig. [1] correspond to the focus-point region. It is worth noticing that, 
unlike the coannihilation or the rapid-annihilation region, the focus-point strip extends to 
remarkably high values of tuq and m\/2- For example, for tan/3 = 50 (lower right panel) 
m ~ 3 TeV (5 TeV) at m 1( / 2 = 1000 GeV (2000 GeV). As a result, the sfermion masses in this 
region are much larger than in the funnel or the coannihilation strip. Hence, for gravitino 
masses up to 1 TeV all the fermion-sfermion decay channels G — > f f are closed, and the 
lifetime is larger than the upper panels. On the other hand, when 1713/2 = 5 TeV the dominant 
decay channels g g, f f and x + W~ are open along the WMAP strips also in the focus-point 
region, resulting in relatively flat lifetime curves as functions of vri\/2- 

We note that for 7713/2 = 5 TeV, the largest value shown, the gravitino lifetime ~ few 
x 10 2 s in all the cases shown in in Fig. [TJ which is comparable with the duration of BBN. 
Lighter gravitinos would decay after BBN is completed. 

As 7711/2 increases for fixed m 3 / 2 , the gaugino masses increase. Therefore, one by one the 
various gravitino decay channels are closed. The last to be closed are the two-body decay 
to Xi 7 an d the three-body decay to light quark pairs Xi QQ- Eventually all the channels are 
closed when m 3 / 2 < m x . For m 3 / 2 = 250 GeV this occurs for m 1( / 2 ~ 580 GeV, as can be 
seen in Fig. [TJ When the dominant channel for hadron production, x\ Z, closes we observe a 
significant decline in the nucleon spectra. This is the reason that in Fig. [3] in all the panels, 
for 7723/2 = 250 GeV the value of B p becomes smaller than 10~ 2 at m\/2 = 400 GeV, and the 
same at mx/2 = 940 GeV when tan/? = 50 and m 3 / 2 = 500 GeV (right panels). Above these 
values of m^, the only channel that produces some hadrons is G — » XiQQ- 

The importance of the channel G — ► x° Z for producing nucleons can be seen in Fig. [2j 
There we plot the quantity e N p (e) for the case of protons for tan/3 = 10 and m 3 / 2 = 
250 GeV (500 GeV) in the left (right) panel. The curves in these figures correspond to the 
various m 1( / 2 we sample along the WMAP coannihilation strip. The reason for the peaks 
in the GeV region is that the protons that are produced by on-shell hadronic decays of 
the Z boson have typical energies of a few GeV. As discussed earlier, for tan/3 = 10 and 
^3/2 = 250 GeV the Xi Z channel closes above m\/2 = 400 GeV. Therefore, we observe two 
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kinds of curves in Fig. [2] (left). These that peak in the GeV region are fed by the XiZ 
channel, whereas these without the peak originate from the three-body channel x? <?<?• F° r 
tan/5 = 10 and m 3 / 2 = 500 GeV, the decay G — > Xi Z * s n °t dosed anywhere along the 
WMAP strip, so all the curves in Fig. [2] (right) exhibit the Z-boson peak. 

There are a few other features in Fig. [3] to be discussed. For my2 ~ 240 GeV in the 
upper panels the lightest chargino can decay on-shell to \i W + . The same happens in 
the lower left (right) panel for my 2 = 520 GeV (m^ = 360 GeV). This causes a sudden 
increase in the number of the produced protons, that is more noticeable for tan/3 = 10 
(upper left panel), since there the relative importance of the decay channel G —>■ xf* W~ — > 
Xi W + W~ is greater. This channel is closed for larger values of rnyz as the chargino becomes 
heavier, resulting in the shoulders we observe for tan/3 = 10 (upper left panel) in the B p 
curves at m 1/2 = 520 GeV (m 3/2 = 500 GeV) and m 1/2 = 800 GeV (m 3/2 = 750 GeV). The 
corresponding features appear also at the same points for tan /3 = 50 in Fig. [3] (upper right 
panel). In the same figure, for m 3 / 2 > 500 GeV, the additional wiggle seen at my 2 ~ 280 
GeV is due to the closing of the G — > xtW~ channel. In the focus-point figures (lower 
panels), at m^j = 400 GeV and m 3 / 2 = 250 GeV the Xi Z channel closes. The same occurs 
at my 2 = 960 GeV (m 1/2 = 1460 GeV) for m 3/2 = 500 GeV (m 3/2 = 750 GeV). After this 
point, B p diminishes. Along the focus-point strip, we also can see in the B p curves the 
effect of the closing of the X12 W~ channels. In particular, the Xi W~ channel closes at 
my 2 = 780 GeV for m 3/2 = 500 GeV, tan/3 = 10, at m 1/2 = 1200 GeV for m 3/2 = 750 GeV, 
tan/3 = 50, and at my 2 = 1640 GeV for m 3/2 = 1000 GeV, tan/3 = 50. For tan/3 = 50, 
the xt W~ channel closes at m 1 / 2 = 440 GeV for m 3 / 2 = 500 GeV, at mi/2 = 760 GeV for 
m 3/2 = 750 GeV, and at mi/ 2 = 1040 GeV for m 3 / 2 = 1000 GeV. These thresholds produce 
distinctive features in the corresponding curves. For very heavy gravitino masses, such as 
the case m 3 / 2 = 5TeV considered here, the dominant two-body decay channels g g, tt, x® Z 
and x + W~ are kinematically available, even in the focus-point region, so no specific features 
are observed, and the nucleon fractions are similar in all the scenarios studied. 

Finally, we note that various other channels close as mi/ 2 increases, without producing 
any significant features in the B p curves. For instance, for tan/3 = 10 (upper left panel in 
Fig. [3]), for m 3 / 2 = 250 GeV, the xV 1 channel closes at mi/ 2 = 340 GeV. The same occurs for 
the channels gg and X\H~ when m 3 / 2 = 500 GeV at mi/ 2 = 220 GeV , and at mi/ 2 = 300 GeV 
for the channel xtW~ . Similarly, for m 3 / 2 = 750 GeV at 77H/ 2 = 260 GeV, the channel xtH~ 
closes, and at mi/ 2 = 340 GeV the channels gg and x~i + H~ close. Just to complete the list 
for this gravitino mass, at mi/ 2 = 400 GeV the Higgs boson channels Xi A H close and at 
Tni/2 = 520 GeV the chargino channel xtW~ closes. 
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5 Constraints on Metastable Particles 



5.1 Generic Constraints on Abundances of Metastable Particles 

Before discussing constraints on specific supersymmetric models with a metastable gravitino, 
we first discuss constraints on the possible abundance (x of a generic metastable particle 
X, as a function of its possible lifetime Tx- We discuss exclusively the constraints due to 
the effects of the electromagnetic and hadronic showers produced in X decays, postponing 
to a later paper discussion of the extra constraints that are imposed on the abundances 
of charged metastable particles by their catalytic effects on light-element abundances due 
to the formation of bound states. Thus, the constraints presented in this Section apply 
exclusively to neutral metastable particles including, but not limited to, the gravitino. Not 
wishing to commit to any specific model, we give results for two typical values of Bh, which 
are applicable also outside the context of specific supersymmetric models. Recall, however, 
that, as noted above, in specific supersymmetric models one may calculate B h and it is not, 
in general, constant. 

The hadronic decays of metastable particles X affect BBN in different ways, depending 
on the stage of BBN in which the non-thermal decay particles interact with the background 
thermal nuclei. This effectively divides the decay effects according to the decaying particle's 
lifetime tx- 

• Decays much before weak freeze-out (tx "C 1 sec) produce showers that are thermalized 
before BBN commences. These decays thus have no impact on light-element abundances, 
and BBN offers no strong constraints on such short-lived decays. 

• Decays during weak freeze-out (1 sec < tx < 100 sec) introduce new interactions that 
may interconvert the neutrons and protons, e.g., via nn + — > pn . These interactions prolong 
the n «-> p equilibrium, and hence delay the freeze-out of the n/p ratio. In this regime the 
effect on the light elements is somewhat similar to the addition of relativistic species, with 
the dominant effect being that on 4 He. However, because the upper limit on 4 He is weak [66], 
this effect does not induce a constraint on our particle properties. 

• Finally, decays following weak freeze-out (tx ^> 100 sec) generate electromagnetic and 
hadronic showers that induce new (photo)nuclear interactions, which in turn may modify 
the light-element abundances established previously by BBN. 

These processes lead to the constraints seen in Figs. SH3 These plot abundance contours 
as a function of pre-decay X abundance (x as in eq. ( fTOl) and lifetime tx, all for decay spec- 
tra corresponding to (rai/2, tan/3) = (300 GeV, 500GeV, 10). In these and subsequent 
figures, the white regions in each panel are those allowed at face value by the light-element 
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abundances reviewed in Section [2TTI Specifically, these are: D/H < 3.2 x 10 5 , 3 He/D < 1.0 
(also shown are dashed lines for 3 He/D < 0.3), Y 4He > 0.240, 6 Li/ 7 Li < 0.05, and 7 Li/H 
< 2.75 x 10~ 10 . The latter constraint comes from 7 Li in globular clusters, and we also de- 
marcate by a dashed line the region within the white area where 7 Li/H < 1.91 x 10~ 10 as 
inferred from field stars. The regions shaded yellow, red, and magenta regions correspond 
to progressively larger deviations from the central values of the abundances as noted on the 
figures. Fig. H] shows the constraints when the effects of both hadronic and electromagnetic 
decays are included. To give a sense of how the various decay modes contribute, Fig. [5] shows 
constraints when hadronic effects are omitted, and only electromagnetic decays are included. 
Figs. M (Fig. Ej) both omit effects of electromagnetic decays, and include only hadronic decays 
to neutrons (protons). 



tan/?=10 m 3/ 2 = 500 GeV 
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Figure 4: Plots of abundance versus lifetime for metastable particles X with lifetimes tx 
between 1 and 10 10 sec, assuming the decay spectra calculated for (m^, ^3/2> tan/?) = 
(300 GeV, 500 GeV, 10), in which case B p m 0.2 and the electromagnetic branching is 
B EM m 3 / 2 = 115 GeV. The X abundance before decay is given by (x = m xnx/n 7 (eq. fTTJj) . 
The white regions in each panel are those allowed at face value by the ranges of the light- 
element abundances reviewed in Section \2.1l whilst the yellow, red and magenta regions 
correspond to progressively larger deviations from the central values of the abundances. 

The basic features of these plots are similar to those found in previous work, but for 
completeness we summarize them here. In all plots, we see that at low (x, all constraints 
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Figure 5: As in Fig. but with electromagnetic decay products only. All hadronic showers 
and resulting interactions with light elements are ignored. 
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Figure 6: As in Fig. [^} but only with decay neutrons. Decay protons and electromagnetic 
particles are ignored. 
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Figure 7: As in Fig. but only with decay protons. Decay neutrons and electromagnetic 
particles are ignored. 

are satisfied except for 7 Li. This is reasonable, as in the limit of Qx — ¥ we recover the 
standard BBN abundances, which agree well with observations except for the 7 Li problem 
which persists. Also, we see a similar behavior in all plots for small tx ^ 10 2 sec. Here, the 
decays occur after weak freezeout but before light element formation occurs, and so for the 
most part the abundances are unaffected; the main exception is perturbations in 4 He due to 
pion interactions delaying n <-> p freezeout. 

For several elements the basic trends at Tx ft 10 2 sec are relatively simple. We begin with 
4 He, which is the only species for which decays always lead to reduced abundances. This is 
physically reasonable, since 4 He is the most abundant complex species, and the non-thermal 
reactions represent sinks but not sources. That is, photoerosion and spallation destroy 4 He, 
but X decays cannot produce it; hence 4 He drops as (x increases. On the other hand, we see 
that decays increase the D/H abundance, which is readily understood physically. Destruction 
of 4 He produces deuteron fragments, some of which are thermalized before interacting and 
thus survive. Because 4 He is so abundant, even if only a small fraction of 4 He is destroyed, the 
resulting deuteron production can be significant. For 6 Li/ 7 Li, the basic trend is also towards 
increasing production with increasing (x- Here, not-yet-thermalized mass-3 fragments can 
interact with ambient 4 He to produce 6 Li via 3 He(a,p) 6 Li and t(a,n) 6 Li. Again, even if 
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only a small fraction of A = 3 nuclei interact this way, this can represent a substantial 6 Li 
abundance compared to observational limits. Moreover, while other secondary processes also 
contribute to 7 Li, the much larger effect is for 6 Li, so that the 6 Li/ 7 Li ratio increases with 
(x- 

By comparing the full constraints (Fig. H]) with those in the electromagnetic-only case of 
Fig. El we see illustrated the well-known result that the electromagnetic decays dominate at 
Tx ?t 10 6 sec, but are ineffective at smaller times [29]. Again for D/H, 7 Li/H, and 3 He/D, 
we see that the hadronic effects from neutrons and protons (Figs. El and d respectively) are 
broadly similar, though the neutron constraints are generally more restrictive out to about 
Tx ~ 10 6 sec. As we will see in detail below, this is the timescale for neutron decay to 
outpace other neutron interactions. 

The cases of 7 Li and 3 He/D are more complicated. Turning to 7 Li in Fig. HI we see that at 
intermediate lifetimes (~ 10 4 sec), decays lead to higher 7 Li/H; this is due to secondary pro- 
duction from unthermalized mass-3 spallation products, 3 He(a,7) 7 Be and t(a, j) 7 Li. How- 
ever, at lifetimes around tx ~ 10 2 — 10 3 sec, there a narrow "valley" emerges in which 7 Li/H 
is reduced, indeed enough to come into agreement with observational limits. By comparing 
Figs. E] and [3, we see that this effect is entirely due to injected neutrons, and is absent when 
only proton decays are included. This suggests that neutrons are destroying 7 Be, and that is 
indeed the case. As noted already by Jedamizk and in subsequent work [21,22], thermalized 
neutrons can destroy 7 Be via 7 Be(n,p) 7 Li and followed by 7 Li(p, a) 4 He. We include this 
effect, but also, as we discuss below (Appendix IA~|) . we allow for further 7 Be destruction 
in which the same reactions are initiated by injected neutrons not yet thermalized. This 
addition slightly enhances the low- 7 Li/H "valley." Notice that the left side of the "valley" 
coincides almost exactly with the constraint from D/H, so that we find D/H > 3.2 x 10~ 5 in 
essentially the entire region where 7 Li/H < 2.75 x 1CT 10 . However, along the left side of the 
"valley" in Fig. HI there may be a marginal solution to the 7 Li problem, to which we return 
later in the context of the CMSSM. 

Finally, turning to 3 He/D at low lifetimes, we also see a reduction in the ratio for large 
(x, with a particularly large effect in the neutron-only case. As with 7 Li destruction, here 
the free neutrons preferentially capture on t and 3 He. This leads to a small 3 He/D ratio as 
(x increases. Specifically, we see that 3 He/D < 0.3 in the narrow strip in Figs. 0] and El where 
the 7 Li may (almost) be solved. 

For a more quantitative understanding we now apply the analytical model developed 
above. 
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5.2 Non-Thermal Abundance Perturbations: Analytical Model 

For the reaction hb — > the rate of decrease of the abundance of the "target" background 
species b is equal and opposite to the rate of production of background species £, namely 

- d t Y b = +d t Y e = Y b T hb ^ e , (27) 

where the reaction rate per target b is given in the thin- and thick-target limits in ( 1251) . The 
net change in these abundances due to this reaction is the time integral of the rate: 

- AY b = +AY e = J Y b T hb ^ £ dt « Y b T hb ^r x . (28) 

In the last expression we take the "exposure time" to be the X lifetime tx = l/r^j this 
is the equivalent of the instantaneous-decay approximation used in the EM case. Since the 
non-thermal spectra, and thus non-thermal reaction rates, scale as Nh ~ Thb-*e ~ Tx, the 
lifetime dependence drops out in the abundance changes. 
For the thin-target case, we have 



AY e 



thin 

bg 



^£ ! Q h (e) ^^de (29) 



Yx%Qh{><*)(?^). (30) 



b_ 

\ n-/,/,_i, J( .| 



where the spectrum-weighted average is 



/ / v L th Qh{e)(T b {e)/a k {e)de 
J eth Qh{e) de 

We see that the abundance change is determined by several parameters. The scaling with Y x 
is intuitively clear-the perturbation should be proportional to the number of non-thermal 
decays per baryon. Another important parameter is the number Qh{> £th) ~ Bh of non- 
thermal h particles per decay above threshold, which is to an excellent approximation the 
total number of h particles per X decay, also an intuitive result. Finally, the cross-section 
factor is an analogue of a branching ratio for the reaction in question relative to all inelastic 
reactions. 

For the thick-target case, we have 

l thick _ v v h g f n /_ \ Qhb-.i/m p _ 



Ay,|JK = wr^ao^* (32) 

= Y x Y? l Q h (> 6,,,) (R)(cu^e/m r ) (33) 
* 6 x 1( n< y xY f (^) (A) (_^_) , (34 ) 
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where we again define weighted averages of the cross section (a) and stopping range (R) 
(with units [g/cm 2 ]) as 

L th Qh(> e) ctm^i (de h /dR)- l de 
{a) ~ f eth Q h (> e) (de h /dR)-ide ' {6b) 
j th Q h (> e) {de h /dR)^de 

{R) = o^) ■ (36) 

Recall that (R) oc Y~± and the fiducial value we choose is appropriate for Y~± ~ 1. We 
see that the abundance perturbations for hadronic decays scale with-and thus constrain-the 
number of pre-decay X particles (gravitinos, in the present study) per baryon (or equivalently 
per photon or per unit entropy). However, to show hadronic as well as electromagnetic results 
on the same plot, it is useful to introduce (x = rnxYxV as in eq. (fTUl) . so that Yx = Cx/v^x- 
We now apply the observed abundance constraints and express them in terms of (x- 

Consider some light-element abundance constraint Yj im . The perturbation saturates this 
constraint when AYe = 5Y° bs = Y^ m — F/ td , where Y^ std is the standard BBN result (see 
e.g. [6]). This abundance will be reached for some value of Yx and thus (. Namely, we have 

<*5L ^ mrfS^/^-V 1 (37) 



KT 11 GeV 




0.5 



{&hb->e/&hk->mel) 



where the fiducial values are appropriate for constraints based on D/H, produced mostly via 
4 He spallation. Here we take the total inelastic interactions hk — ► inel to be dominated by 
interactions with background 4 He, so that Y^ e /Y^ g = 1; this is the case for nucleons up to 
e ~ 10 GeV, and beyond this the smaller NN inelastic cross sections change the result by 
no more than about a factor of 2. 

We compare our analytical estimates to numerical results for the production of deuterium 
obtained using our codes, as shown in Figs.HHTl We see that for short lifetimes, our prediction 
from eq. ( [371) agrees with the D/H abundance constraints for Yd = 3.2 x 10~ 5 to within a 
factor of 3. In particular, our thin-target estimate is a good description of the neutron-only 
results (Fig. [6]). This was to be expected as neutron propagation is, for these X lifetimes, 
well within the thin-target regime. 
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For the thick-target case, putting AYt = 5Y° hs gives 
C mck ~ ^° bs ^-i rn p /(a hb ^ e ) 

= 10~ 9 GeV (39) 

x 




where again we use fiducial values that are appropriate for constraints based on D/H, pro- 
duced mostly via 4 He spallation and Y e ± = 1. We see that the thick-target limits are about 
a factor ~ 100 weaker than the thin-target limits. This is reasonable, as the thick-target 
limit is less efficient in light-element production. Again comparing Figs. HE we see that the 
proton-only results (Fig. [7j) lie between the thin- and thick-target estimates. This is also 
as expected: for the X lifetimes of interest, the proton propagation is partially within the 
thick-target regime. 

Given our success in understanding of D/H, it is worth extending the analysis to 6 Li, 
which should be due to secondary production of non-thermalized mass-3 nuclides via 3 He(a, p) 6 Li 
and t(a,n) 6 Li. If this is the case, then the 6 Li abundance should be given by the product 
of the non-thermal mass-3 abundance AY 3 and the probability that the mass-3 nuclide in- 
teracts before it stops. But for the thin-target case, the probability is small, and is just 
given by the product of the mass-3 interaction rate T 3a = Y^n^a^aV times its stopping time 
r 3,sto P = (/ b" 1 de) ps (m p n B v / R?,)" 1 , with R 3 the mass stopping grammage. Thus we have 

AY(«Li) = AY 3 Y^R ^ - N>U \ (40) 

Trip 

with Y^ g ~ 0.1. Using the above thin-target expression for AY d (eq. [291 . and taking 
( 3 He/D) spa ii ~ 1 for spallation production, this becomes 

AYTLi) „ S^bJ^-) ^ ^EL^L . (41) 

7]m X \ Vhk^mcl I fn p 

Assuming the 7 Li perturbation is small, so that Y(jLi) ps l A ( 7 Li)BBN,std ~ 3 x 10 -10 , we solve 
for (x in terms of the observed limit on 6 Li/ 7 Li = y( 6 Li)/y( 7 Li): 

(St M * ^( 6 W 7 Li) obs ^^(^)- 1 R ^ A :i N , u) (-) 
~ 2 - 5Xl ° GeV [ 0.05 V ^ J L(3Aa^N°U) ' (43) 
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where the fiducial values not shown are the same as in eq. (1371) . This result within a factor 
of 3 the strongest 6 Li/ 7 Li constraints around tx ~ 10 4 sec, where the neutrons have time to 
interact before they decay. 

The agreement between our order-of-magnitude estimates and full numerical results 
serves as a strong check on our analysis and gives us confidence in both approaches. 

5.3 Constraints on the Abundance of a Massive Gravitino 

We now apply the above analysis to the specific class of supersymmetric scenarios with a 
metastable massive gravitino and a neutralino LSP \. 

We recall that neutralino LSP scenarios are characterized by narrow strips along which 
the relic x density lies within the range 0.0975 < ficDM^ 2 < 0.1223 that is favored by 
WMAP and other observations [1]. Since this range is relatively narrow (namely only a 
few %), within any specific supersymmetric scenario it determines some combination of the 
model parameters also to within O(few) %. A more detailed discussion is given in [47], 
where some explicit parameterizations of such WMAP strips are given. In the examples 
given there, the value of rriQ is tightly determined in CMSSM scenarios as a function of 
for fixed values of tan/3 and Aq. Many properties of these supersymmetric models change 
little as one varies mo across such a narrow WMAP strip. Specifically, the branching ratios 
for massive gravitino decay and hence vary little across a strip, and one may usefully 
represent the cosmological constraints on such CMSSM scenarios as functions of m\ii alone, 
using a representative value of Bh that is calculated as a function of mi/2- 

One representative example is shown in Fig. [HJ the first five panels, with shadings, display 
the effects of the decays of a gravitino with a mass 777.3/2 = 250 GeV on the different light- 
element abundances (D/H, 3 He/D, 4 He, 6 Li/ 7 Li and 7 Li/H) as functions of mi/2 along the 
WMAP strip in the coannihilation region for a CMSSM scenario with tan/3 = 10, Aq = 0. 
As in Fig. H] - [TJ the white regions in each panel are those allowed at face value by the 
ranges of the light-element abundances reviewed in Section 12. 1\ and the yellow, red, and 
magenta regions correspond to progressively larger deviations from the central values of 
the abundances. The final panel (bottom right), shows the strongest constraints from each 
abundance, as labeled, and from these one infers the strongest overall constraint, shown by 
the thick black curve. 

We first note that all constraints weaken abruptly as mi/ 2 — > 600 GeV. This corresponds 
to the limiting case when m x ~ QA2mi/ 2 — > 7773/2 = 250 GeV. In this limit, the energies 
in the EM and HD showers vanish, and the effects of gravitino decay disappear. Note 
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Figure 8: The effects of the decays of a gravitino with a mass m.3/2 = 250 GeV on the different 
light- element abundances (D/H, 3 He/D, 4 He, 7 Li/H and 6 Li/ r Li y ) as a function ofrrti/2 along 
the WMAP coannihilation strip for a CMSSM scenario with tan/? = 10,v4 = 0. As in 
Figs. [7} the white regions in each panel are those allowed at face value by the light-element 
abundances reviewed in Section l2J\ and the yellow, red, and magenta regions correspond to 
progressively larger deviations from the central values of the abundances. 



that mi/ 2 = 250 GeV is the smallest value we consider; for larger values this effect is less 
important. 

Consider now the D/H constraints in the region m\/2 ^ 600 GeV in the top left panel of 
Fig. [SJ We see that the lowest contour, at D/H = 3.2 x 10~ 5 , lies at about ( 3 / 2 ~ 10~ n GeV 
between mi/ 2 = 180 GeV and just under 600 GevQ. To understand this behavior, we recall 
from Fig. [T]that for mi/ 2 = 180 — 580 GeV, the gravitino lifetime grows from about 10 7 sec 
to 10 10 sec, and in Fig. Owe see that the nucleon branching ratios B p w B n ~ 6 x 10~ 2 for 
m i/2 400 GeV. Combining these with the generic lifetime dependence in Fig. HJ we see 
that for B p w B n ~ 2 x 10 _1 , the lowest D/H contour for r > 10 7 sec is roughly constant 
at (3/2 ~ 3 x 10~ 12 GeV. In our case, these constraints weaken due to the lower branching 
ratios by a factor ~ 3, yielding the value (3/2 ~ 10~ n GeV seen in Fig. [8] The higher D/H 
contours in Fig. [8] can be understood in a similar manner. 

As another example, consider the 7 Li/H panel in the bottom middle panel of Fig. [HJ 

1 We recall that the gravitino becomes the LSP for m 1 / 2 <; 600 GeV, in which case a different analysis is 
necessary. 
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There we see, for m.1/2 < 600 GeV, an "island" of the lowest constraint where £3/2 ~ 
10~ 10 — 10 -8 GeV. As (3/2 drops, 7 Li/H then grows, peaking to highest level where C3/2 ~ 
3 x 10~ 9 GeV. Then as C3/2 continues to decrease, 7 Li/H continues drops down to its standard 
BBN level. From the generic plot Fig. HI we see that for r > 10 7 sec and with B p ~ 0.2, the 
7 Li/H constraint forms an island around C3/2 ~ 10~ 9 GeV, with a width that grows with r. 
Scaling down to the lower branching, we again can understand the behavior in Fig. [HJ 

The fact that this qualitative analysis works so well indicates that the generic results 
are a good approximation to the full ones, modulo changes in branching ratios. However, 
the generic results are for a particular decay spectrum whereas, as we have seen, there are 
significant variations in the decay spectrum. Thus we infer that our results are not strongly 
dependent on the detailed shapes of the decay spectra beyond the sensitivity to the branching 
ratios. This also makes sense in terms of our analytic approximations, which suggest that the 
decay spectra enter principally via their integral properties and particularly their branching 
ratios. 

This understanding of the non-thermal particle effects leading to our constraints gives 
confidence in our results, and thus to their implications for supersymmetry. Namely, again 
in the case of Fig. [HJ where tan/? = 10 and m 3 / 2 = 250 GeV, the fact that the observed 
abundances generally agree with the standard BBN calculations implies that the allowed 
(white) regions are generally at low gravitino abundance. The exception is the bottom 
middle panel, where the observational discrepancy with the standard BBN calculation of 
the 7 Li abundance is reflected in the fact that the yellow region extends from a vanishing 
gravitino abundance up to quite large values, and the "preferred" white region appears only 
when £3/2 > 10~ 9 GeV and m.1/2 is not too large B It is clear that this "preferred" region for 
7 Li/H is incompatible with the allowed regions for the other light-element abundances shown 
in the other panels. Thus, there is no value of m.1/2 in this particular CMSSM scenario that 
solves the 7 Li/H problem. As already mentioned, we disregard the 7 Li/H problem in the 
compilation of constraints shown in the bottom right panel of Fig. [HJ and in the following 
discussion. 

The weakening of the constraints as m.1/2 — » 600 GeV implies that a large abundance is 
allowed for m 3 / 2 = 250 GeV if m-i/2 ~ 600 GeV, apart from the issue of the 7 Li abundance. 
We recall that if m.1/2 > 600 GeV with fixed m.3/2 = 250 GeV, the gravitino becomes the 
LSP and the lightest neutralino becomes the NLSP. In this case, the constraint due to the 
cosmological relic density should be applied to the gravitino, the WMAP strip is no longer 

2 We also note the existence of a more disfavored (red) region in the bottom middle panel showing the 
7 Li/H ratio, appearing when £ 3 / 2 > GeV and m 1 / 2 ~ 500 GeV. 
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relevant, and a different analysis would be required. 

In the final, summary panel of Fig. [8] (bottom right), we see explicitly that the 7 Li 
constraint (light blue) is incompatible with the other constraints due to 4 He (green), D/H 
(magenta), 3 He/D (red) and 6 Li/ 7 Li (dark blue). The weakest of these constraints is that due 
to 4 He, the two strongest constraints at smaller and larger m^, respectively, are those due 
to 6 Li/ 7 Li and 3 He/D, and the combined constraint is shown in black. We see that it hovers 
in the range C3/2 = ^3/2^3/2/^7 ~ 3 x 1CT 12 to 1CT 12 GeV, corresponding to n 3 / 2 /n 1 ~ 10~ 14 , 
except in the limit as — > 600 GeV. 

Figs. [HI [TOl [HJ and [12] show the correspondingly coloured regions of varying discomfort 
for the different light-element abundances for the larger values of 7723/2 = 500 GeV, 750 GeV, 
1000 GeV and 5000 GeV, respectively. The constraints are shown for the full length of the 
WMAP strip up to m 1 / 2 ~ 900 GeV, along which m x < m 3 / 2 for all the displayed values of 
7723/2- For m 3 /2 < 1000 GeV, we see that the D/H constraint (top left panel in each plot) 
and the 4 He constraint (top right in each plot) are relatively stable, reflecting the rough 
constancy of the D/H results for long lifetimes (see Fig. HJ), and the rough constancies of the 
branching ratios (see Fig. [3]). For m 3 /2 = 5000 GeV, the short lifetimes (see Fig. [IJ) severely 
weaken the constraints for all elements (cf. Fig. HJ), with D/H surviving as the strongest for 
all mi/2- We discuss later the case of large m 3 / 2 . 
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Figure 9: As for Fig. 0, with m 3 / 2 = 500 GeV but unchanged values for the CMS SM param- 
eters. 
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Figure 10: As for Fig. 0, with 777,3/2 = 750 GeV but unchanged values for the CMSSM 
parameters. 
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Figure 11: As for Fig. 
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with 7773/2 = 1000 GeV but unchanged values for the CMSSM 
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Figure 12: As for Fig. 0, with m 3 / 2 = 5000 GeV but unchanged values for the CMSSM 
parameters. 

In the case of 7 Li/H (bottom middle panels), the "preferred" white region changes posi- 
tion as increases, moving to larger 777i/ 2 for 7713/2 = 500, 750 GeV, but reverting to low 
77ii/2 a ^ high 777.3/2- The mass ranges are those which have lifetimes near ~ 3 x 10 6 sec 
(see Fig. 0]), whereas the 7 Li constraint is strongest for long lifetimes. However, there is 
no overlap between the white regions in this and the other panels, implying that the 7 Li 
problem cannot be solved for any value of m 3 / 2 for the particular values of tan/3 and A 
chosen here. As stated previously, we do not include the 7 Li constraint in our compilation. 
As 7 Li is problematic in standard BBN, the inclusion of this constraint would give the false 
impression that nearly every supersymmetric model with a decaying gravitino is excluded. 
We are therefore implicitly assuming that there is another solution for the 7 Li problem, e.g., 
due to observational or astrophysical uncertainties as discussed in Section |2~T1 In contrast, 
the constraints from the other light elements perturb the previous concordance of BBN with 
respect to those elements. 

Finally, we note that for m 3 / 2 = 500 to 1000 GeV, the 6 Li/ 7 Li constraint is significantly 
stronger than for m 3 / 2 = 250 GeV, becoming progressively stricter. This reflects the stronger 
limits arising when the gravitino lifetime is shorter (see Fig. H]). In all cases the 6 Li/ 7 Li ratio 
provides the most restrictive limit on £ 3 / 2 , and strengthens by a factor ~ 10 as 7771/2 increases 
from 250 to 1000 GeV. Specifically, we find C3/2 ^ 10" 12 - IO" 11 GeV for m 3/2 = 500 GeV, 
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Cs/2 < 3 x 1(T 13 - 3 x 1(T 12 GeV for m 3/2 = 750 GeV, and Cs/2 < 
m 3/2 = 1000 GeV. 

In the case of Fig. [12], we see that the most significant upper limit on the gravitino 
abundance is that from D/H. Comparing this with the lower limit on the gravitino abundance 
coming from 7 Li/H, we see that they are marginally compatible over essentially the full 
range of displayed. However, this conclusion is crucially dependent on the precise 
implementations of the D/H and 7 Li/H constraints: for example, if the upper limit on the 
7 Li/H abundance is strengthened to 1.91 x 10~ 10 , as suggested by field stars and indicated 
by the dashed line in the middle lower panel of Fig. [12], compatibility becomes more difficult. 

Figs. [TBHTTI are the same as Figs. IBHT21 but for tan (3 = 50. Similar trends emerge as for 
the tan /3 — 10 results, with some differences of detail due to the more rapid rise in lifetimes 
with mx/2 (Fig. [[]). The range in mi/ 2 is extended to almost 1.9 TeV and includes the rapid 
annihilation funnel in addition to the coannihilation region. We again see that no regions 
allow for a solution to the Li problem while simultaneously satisfying the other light-element 
constraints. The constraints on £ 3 / 2 are again dominated by the 6 Li/ 7 Li ratio, and are weak 
for m 3 / 2 = 5000 GeV. As in the tan/5 = 10 case, the 7 Li/H constraint is incompatible with 
the others for m 3 / 2 < 1000 GeV, but may be marginally compatible for m 3 / 2 = 5000 GeV. 
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Figure 13: As for Fig. [Sj, with m 3 / 2 = 250 GeV but tan/3 = 50 and unchanged values for the 
CMSSM parameters. 

Figs. [T8H221 present a similar analysis to that in Figs. |8HT2"| but for the focus-point region 
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Figure 14: ^4s for Fig. 0, ^3/2 = 500 GeV, but tan/5 = 50 and unchanged values for 
the CMSSM parameters. 



10- 7 
10" 8 

" 10 " 9 
£ 10"'° 

j« 10"" 

1Q -.2 

10"' 3 

10" 7 
10" 8 

^ 10" 9 

§ io-'° 

- 10"" 

to 

10" 12 
10"' 3 

io-' 4 



tan/S=50 m 3/ 2 = 750 GeV 

















■ 








0.240 V 


■ 3.2 x 10/ ' 








D/H 


3 He/D 


"He 






\ N'LI/H \ 
3 He/o\ 








2.75 xlO' 11 10 * 10 
















\i/ 7 Li 

... 1 .... 1 .... 1 ■ 




7 Li/H 




Ai/'L! 
oil 











500 10001500 
m I/2 [GeV] 



500 10001500 
m 1/2 [GeV] 



500 100015002000 
m 1/2 [GeV] 



Figure 15: ^4s /or Fzg. 0, wtt ^3/2 = 750 GeV, but tan/3 = 50 and unchanged values for 
the CMSSM parameters. 
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Figure 16: As for Fig. 0, with (top left) 777.3/2 = 1000 GeV, but tan/? = 50 and unchanged 
values for the CMSSM parameters. 
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Figure 17: ^4s /or Fig. 0, to^ 7773/2 = 5000 GeV, but tan/3 = 50 and unchanged values for 
the CMSSM parameters. 
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of the CMSSM with tan (3 = 10. The results are also largely similar to those for the previous 
cases. The same is true for the focus-point region of the CMSSM with tan/? = 50, shown in 
Fig. [2314271 As in the previous cases, the 7 Li/H constraint is incompatible with all the other 
constraints, except possibly for m 3 / 2 = 5000 GeV. 
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Figure 18: As for Fig. [21 with ^3/2 = 250 GeV, but tan/5 = 10 and CMSSM parameters 
appropriate for the WMAP strip in the focus-point region. 




Figure 19: As for Fig. [21 with 1113/2 = 500 GeV, but tan/? = 10 and CMSSM parameters 
appropriate for the WMAP strip in the focus-point region. 
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Figure 20: As for Fig. 0, with m 3 / 2 = 750 GeV, but tan/3 = 10 and CMSSM parameters 
appropriate for the WMAP strip in the focus-point region. 
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Figure 21: v4s /or i<%. 0, with m 3 / 2 = 1000 (7eV, but tan/5 = 10 and CMSSM parameters 
appropriate for the WMAP strip in the focus-point region. 
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Figure 22: As for Fig. wift m 3/2 = 5000 GeF ; 6u< tan/3 
appropriate for the WMAP strip in the focus-point region. 
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Figure 23: ^4s /or [21 ^3/2 = 250 GeV, but tan/3 = 50 and CMSSM parameters 
appropriate for the WMAP strip in the focus-point region. 
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Focus Point ton/? = 50 m 3/2= 500 GeV 
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Figure 24: As for Fig. 0, with 7713/2 = 500 GeV, but tan (3 = 50 and CMSSM parameters 
appropriate for the WMAP strip in the focus-point region. 
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Figure 25: ,4s /or Fig. 0, wt/i 7773/2 = 750 GeV, but tan/3 = 50 and CMSSM parameters 
appropriate for the WMAP strip in the focus-point region. 
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Figure 26: As for Fig. [21 with m.3/2 = 1000 GeV, but tan (5 = 50 and CMS SM parameters 
appropriate for the WMAP strip in the focus-point region. 
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Figure 27: v4s for Fig. [21 wift m 3 / 2 = 5000 GeV, but tan/? = 50 and CMSSM parameters 
appropriate for the WMAP strip in the focus-point region. 
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5.4 Varying m 3 / 2 and the Li Problem 

The previous plots displayed the gravitino constraints along WMAP strips in the (mu2, m ) 
planes for a few discrete choices of m 3 / 2 . We now display some results as continuous functions 
of for a few WMAP-compatible points with certain discrete values of m^. Fig. [281 
shows our first choice, = 400 GeV and tan/3 = 10, in which case the WMAP point 

is essentially benchmark point C defined in [47]. As expected from the previous figures, we 
see that the other light-element constraints are incompatible with the 7 Li/H constraint for 
low m 3 / 2 . Disregarding the 7 Li problem, we find, e.g., an upper limit on £ 3 / 2 ~ 10 -13 when 
m 3 / 2 ~ 1.4 TeV. However, marginal compatibility with the cosmological 7 Li abundance is 
approached for m 3 / 2 ?t 3 TeV with £ 3 / 2 > 10~ n GeV, as was to be expected from Fig. [121 This 
is a realization within the CMSSM of the marginal compatibility between the 7 Li abundance 
and the D/H ratio that was noted earlier in the context of Fig. HJ 
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Figure 28: The effects of the decays of a gravitino with variable mass m 3 / 2 on the different 
light- element abundances for a specific point (benchmark C) with 7711/2 = 400 GeV on the 
WMAP coannihilation strip for a CMSSM scenario with tan/3 = 10, Aq — 0. As in previous 
figures, the white regions in each panel are those allowed at face value by the light-element 
abundances reviewed in Section \2.1\ and the yellow, red, and magenta regions correspond to 
progressively larger deviations from the central values of the abundances. We see marginal 
compatibility between the 7 Li constraint (light blue) and the other constraints for m 3 / 2 > 
3 TeV. 
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Figs [291 [30] and [31] display similar features for other choices of WMAP-compatible bench- 
mark points [47]. In the case of Fig. [29J we choose benchmark point L with m\i 2 = 460 GeV 
and tan/3 = 50, which is also in a coannihilation strip. We see very similar features to Fig. [28l 
e.g., an upper limit on £ 3 / 2 ~ 2 x 10 -13 GeV when m 3 / 2 ~ 1.6 TeV if the 7 Li constraint is disre- 
garded, and a marginal solution of the 7 Li problem for m 3 / 2 > 3 TeV with £ 3 / 2 > 10~ n GeV. 
Fig. [30] is based on benchmark point M with m^j = 1840 GeV and tan/3 = 50, which 
is in a rapid-annihilation funnel. The features seen in the previous figures shift to higher 
m 3/2, e.g., we see an upper limit on ( 3 / 2 ~ 10~ 12 GeV when m 3 / 2 ~ 2.4 TeV if the 7 Li con- 
straint is disregarded, and a marginal solution of the 7 Li problem for m 3 / 2 ^ 4 TeV with 
C3/2 > 5 x 10~ n GeV. Finally, Fig. [31] is based on benchmark point E with mi/ 2 = 300 GeV 
and tan/3 = 10, which is in a focus-point region. We see in general results that are very 
similar to those for benchmark point C shown in Fig. [28] 
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Figure 29: As for Fig. [28\ for CMSSM benchmark point L with = 460 GeV and 

tan (3 = 50, on the WMAP strip in the coannihilation region. 

The stability of the marginal "solution" of the 7 Li problem reflects the fact that the 
gravitino lifetime ~ 10 2 — 10 3 s and the number of nucleons per decay are relatively stable for 
heavy gravitino masses as mi/ 2 and tan (3 are varied: see Figs [1] and [3] The underlying physics 
of this marginal "solution", as discussed in [21,22], is that for this range of gravitino lifetimes, 
thermalized neutrons can destroy 7 Be via 7 Be(n,p) 7 Li, following which 7 Li is destroyed by 
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Point M tonj8 = 50 m, /2 =1840 GeV 
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Figure 30: As for Fig. H3, for CMSSM benchmark point M with m^j = 1840 GeV and 
tan (3 = 50, on the WMAP strip in the rapid- annihilation funnel region. 
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Figure 31: As for Fig. [28\ for CMSSM benchmark point E with m\ii = 300 GeV and 
tan (3 = 10, on the WMAP strip in the focus-point region. 
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the 7 Li(p, a) 4 He reaction. Here, as discussed in Appendix we also allow for supplementary 
7 Be destruction by non-thermalized neutrons. 

6 Conclusions 

The first objective of this paper has been to document a new suite of codes for the de- 
velopment of non-thermal particle showers in the early universe. These codes treat both 
electromagnetic and hadronic components of showers, and can be applied to the decays of 
unstable particles within a wide range of possible lifetimes, that may decay either during or 
subsequent to standard BBN. We have also shown how this suite of codes may be used to 
analyze the possible effects of generic particles with typical decay modes. 

As a specific application of this suite of codes, we have re-examined the cosmological con- 
straints on unstable gravitinos in /^-conserving CMSSM scenarios with a neutralino LSP aris- 
ing from their effects on the light-element abundances. The first step in this re-examination 
was the recalculation of some two-body gravitino decay modes and the calculation of some 
important three-body decay modes. The next step was the simulation of the electromagnetic 
and hadronic products of these decays using PYTHIA. We have then used our new suite of 
codes to simulate the interactions of these decay products with the cosmological plasma, 
including both energy losses and interactions with the nuclei of light elements that affect 
their abundances. Generally, we find that details of the decay spectra are less important 
than the overall fractions of baryons produced in the gravitino decays. These fractions are 
model-dependent, and vary considerably with the gravitino mass and the CMSSM parame- 
ters. In addition to the baryon fractions, the constraints on the gravitino abundance inferred 
from the observed light-element abundances depend sensitively on the gravitino lifetime, and 
are generally weaker in scenarios with a shorter-lived gravitino. 

For m3/2 < 3 TeV, none of the CMSSM scenarios we study solves the cosmological 7 Li 
problem. Accordingly, one must either reject these CMSSM scenarios, or re-evaluate the 
observational data on 7 Li, or postulate some other mechanism for bringing the 7 Li abun- 
dance into line with standard BBN predictions, or find some way to modify these predic- 
tions. Setting the 7 Li problem aside, we find that the strongest constraints are usually 
those imposed by the 6 Li/ 7 Li ratio, with weaker limits coming from the D/H, 3 He/D and 
4 He constraints. We have evaluated these constraints along WMAP strips in the (rn.1/2, mo) 
plane where the relic neutralino density falls within the expected range for cold dark mat- 
ter, for m 3 /2 = 250, 500, 750, 1000 and 5000 GeV, A = and two representative choices 
tan/3 = 10,50. There are two strips for each choice of tan/3: one in the coannihilation and 
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funnel region, and one in the focus-point region. 

We find that the upper limits on ( 3 / 2 = ^3/2^3/2/^7 are quite similar along these WMAP 
strips, that they weaken as mi/ 2 increases for fixed m 3 / 2 (particularly for tan/5 = 50), that 
they strengthen by about an order of magnitude as m 3 / 2 increases from 250 GeV towards 
1000 GeV, and that the constraints are significantly weaker for m 3 / 2 = 5000 GeV. The 
constraints also weaken considerably as the neutralino mass approaches m 3 / 2 from below, as 
the phase space for gravitino decay disappears. For larger values of m^, the neutralino is 
no longer the LSP, and the WMAP strips are inapplicable. 

Extending these studies to larger m 3 / 2 , we find that the D/H and 7 Li/H constraints 
may become marginally compatible for m 3 / 2 > 3 TeV and a very narrow range of gravitino 
abundance that increases with m 3 / 2 . In this region, both thermal and non-thermal neu- 
trons destroy 7 Be via the 7 Be(n,p) 7 Li reaction, which is followed by 7 Li destroyed via the 
7 Li(jo, a) 4 He reaction [21,22] for a finely-tuned range of lifetimes around t x ~ 10 3 sec. We 
use here the recently updated 3 He(o;,7) 7 Be rate [99] which makes the 7 Li problem worse, 
since it requires 7 Be destruction by a larger factor. We find that the competing upward 
perturbations to D/H and 3 He/D drive these species away from their observed levels in the 
regimes where the 7 Li problem is solved. This highlights the continued importance of the 
nuclear reaction data, particularly at energies beyond the < 1 MeV range important for 
standard, thermal BBN reactions. Additional data extending from ~ 3 MeV to higher en- 
ergies ~ 30 MeV are urgently needed, particularly for reactions such as 3 He(o;, 7) 7 Be, and 
for 7 Be + n — > anything. Further studies of both the light-element abundances as well as 
the relevant nuclear reactions is required in order to verify whether this possible "solution" 
to the 7 Li problem is viable and, if so, how finely one must tune the gravitino abundance 
within the CMSSM. 

The extension of this analysis to the case of a gravitino LSP depends on the nature of the 
NLSP. If the NLSP is neutral, much of the suite of codes developed here would be directly 
applicable. On the other hand, if the NLSP is charged, e.g., the lighter stau or some other 
slepton, it may form electromagnetic bound states before decaying, and catalyze nuclear 
interactions that are not included in the codes developed here. We will return to this issue 
in a future publication. 
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Appendices 

A Propagation of Non-Thermalized Particles in the 
Cosmic Plasma 

We wish to find the non-thermal spectra Nh which are solutions of the propagation equation 
f|T6|) . As discussed in Section [372| the equation is self- regulating, with N h adjusting itself so 
that sources come into equilibrium with sinks, and hence that d t Nh = 0. In this case, the 
equation we wish to solve is 

J h (e) - T h (e)N h (e) - d e [b h (e)N h (e)] = 0. (A-l) 

The sink term has two contributions, due to elastic and inelastic scattering: 

T h(e) = ^2n b [av] hb ^ h (e) + ^ n b [av) hb ^ h >(e). (A-2) 

b h'^h,b 

The source term has three contributions, due to direct injection, elastic down-scattering and 
inelastic down-scattering, respectively, given by: 

Jh\ e ) — Jx,h( e ) + ^elastic./i ( e ) + ^inelastic,h (e) 

= T x Q h {e) + 2^n b J N h — de 

ST [°° at d[av} h ^ h {e',e) 

+ 2^ n *> / Nh ' j e " e ' ( A_3 ) 

h'^h,b ^ e <( 

where gives the number of hadron of type h per X decay as discussed in Section [31 so 

that Jx,h = F x Qh gives the h injection rate per X particle. The second term in eq. (1A-3I) 
gives the rate at which h particles at energy e are produced by elastic scatterings of other 
h particles, integrated over initial energies e'. The third term in eq. (IA-31) gives the rate at 
which h particles at energy e are produced via inelastic processes h'b — > h, integrated over 
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initial energies e' . Note that the elastic and inelastic sources in eq. (1A-3I) are both integrals 
over proton and neutron spectra N p and N n . Thus, in the cascade equation eq. (TIE]) , there is 
one factor of a hadronic spectrum Nh in every term except the decay injection term. Thus, 
on dividing the cascade equation by T x , one infers an equation for N h = N h /T x which has 
the same energy losses and scattering sources and sinks, but whose decay injection term 
is Qh(e), which is independent of Tx- In other words, we find that the physical spectrum 
depends in the following way on the decay properties: N h (e) = r x N h (e), i-e., the scaling 
separates A^(e) into two factors: the X decay rate, and the branching and spectral shapes 
of the h daughters. 

We now look at each of the source terms in detail. 

A.l Elastic- and Inelastic-Scattering Terms 

Hadronic processes are described in ( 1A-2I) and ( 1A-3I) by the differential cross section da^b-*^ e?) / dej 
for hadron hi with energy to produce hj with energy €j. The corresponding total cross 
section for the process is 



We find it useful to split the collisional loss terms into the rates for inelastic collisions with 
background nuclei: X h ncl = J2b n b a hb l (. e ) an d elastic collisions: A^ 1 = ^& n 6 <J ^-*^( e )' where h 
loses energy. Similarly, we have divided the source terms into a double sum over all inelastic 
processes j s hower^b g - > i and a sum over elastic scatterings. 

The kinematics of the final state particles are treated as in [22]. 

A. 2 Energy Losses 

A decay particle loses energy not only loses discontinuously (in time and in energy space) 
through two-body scatterings, but also continuously (in time and energy space) through 
collective electromagnetic interactions with the background plasma. These are the analogs 
of the usual ionization energy losses of particles in laboratory matter, but modified to reflect 
the fully-ionized state of the thermal plasma in the early universe. The collective interactions 
amount to a sum of soft, long-range two-body interactions with the surrounding particles, 
which occur out to some screening radius. As a result, the energy- loss rates are linearly 
proportional to the density of the background species in question, and to the electromagnetic 
coupling for the interaction in question. 




(A-4) 
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Table 2: Nuclear reactions used in this analysis. 



Reactions used in non-thermal propagation 



Reaction 



Reference 



NN -> NN (elastic) 
NN -> AW (total) 
p 4 He — > p 4 He (elastic) 
n 4 He — > p 4 He (elastic) 

— * (inelastic) Meyer 
p 4 He — > (inelastic) 
p 4 He -> d 3 He 
p 4 He — > np 3 He 
p 4 He — * ddp 
p 4 He — >■ cinpp 
p 4 He — > nnppp 
p 4 He -> iV 4 He7r 
n 4 He — >■ (inelastic) 
n 4 He — > npt 
n 4 He — > ddn 
n 4 He — ► cinnp 
n 4 He — > nnppp 
n 4 He -> N 4 Reir 



Kawasaki et al. [22] 
Kawasaki et al. [22] 

Meyer [100] 

Meyer [100] 
[100]; Ando, Cyburt, Hong, and Hyun [101] 

Meyer [100] 

Meyer [100] 

Meyer [100] 

Meyer [100] 

Meyer [100] 

Meyer [100] 

Meyer [100] 

Meyer [100] 

Meyer [100] 

Meyer [100] 

Meyer [100] 

Meyer [100] 

Meyer [100] 



Additional reactions with background nuclides 



Reaction 



Reference 



pn — > d'j 
pd — ■> (inelastic) 
pd — > 3 He7 
pt — > n 3 He 
pt — > (inelastic) 
p 3 He — > (inelastic) 

p 6 Li -> 7 Be7 
p 7 Li -> 4 He 4 He7 
p 7 Be -> 8 B7 
d 4 He -> 6 Li 7 
t 4 He -> 6 Lin 
t 4 He -> 7 Li 7 
3 He 4 He -> 6 Lip 
3 He 4 He -> 7 Be7 
nc? — > (inelastic) 

nrf — > t'j 
nt — > (inelastic) 
n 3 He — > (inelastic) 
n 6 Li -> 7 Li 7 



Ando, Cyburt, Hong, and Hyun [101] 
Kawasaki et al. [22] 
Cyburt et al. [20] 
Cyburt [5] 
Kawasaki et al. [22] 
Kawasaki et al. [22] 
Cyburt et al. [20] 
Cyburt [5] 
this work 
Mohr [102] 
Cyburt et al. [20] 

Cyburt [5] 
Cyburt et al. [20] 
Cyburt and Davids [99] 
Kawasaki et al. [22] 

Cyburt et al. [20] 
Kawasaki et al. [22] 
Kawasaki et al. [22] 
Cyburt et al. [20] 
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Thus, the energy-loss rate bh{e) = —(deh/dt) background of hadron species h sums over the 
non-thermal particle's electromagnetic interactions with the background plasma: 



where the terms account for interactions with background electrons, positrons and photons, 
as well as pair-production. Before going into the details of the energy dependences of these 
terms, some general effects are important to note. 

The dominant electromagnetic losses for both protons and neutrons are due to interac- 
tion with the background electrons and positrons, which are the most mobile and abundant 
charged species. Consequently, the key loss rates are proportional to the total number density 
n e ,tot = n e -+n e + of background electrons and positrons, itself set by the plasma temperature 
and the excess n e ^ et = n e - —n e + = {Z)-qtib ~ (proton + 2l4 He )nB of electrons over positrons 
required by charge balance with background protons and nuclei. As the universe under- 
goes pair annihilation, the electron/positron abundance and thus the electromagnetic losses 
change drastically, and pass through three regimes [103]. At T > m e , pairs are relativistic 
and have numbers comparable to photons: n e> tot ~ 3n 7 /2 ^> n e ,net- At m e > T > m e /25, 
electrons are nonrelativistic but pairs remain much more abundant than baryons, though 
less than photons, with n ejto t ~ \J n l,nct + (2^ e +) 2 and n e + = 2(m e T/2n) 3 ^ 2 e~ me ^ T . At 
T < m e /25, pairs are nonrelativistic, the positron density is vanishingly small, and the 
electron density is set by charge balance: n e tot = n e net . 

Non-thermal ions suffer all of the losses represented by the terms of eq. (IA-50 . while 
neutrons see only the term that accounts for interactions via their magnetic moment. A 
consequence of this distinction is that charged decay hadrons have a qualitatively different 
energy-loss behavior from that of neutrons; we thus discuss the two cases separately. 

For non-thermal ions, the dominant energy losses are due to Coulomb interactions with 
the background electrons and positrons. Such losses of fast ions in plasmas have a large 
literature of dedicated study [104-106] in addition to calculations made in support of early- 
universe applications [11,27]. Following the general approach of [105], we write Coulomb 
losses for an ion with energy 7M and charge Z in the form 




(A-5) 




(A-6) 



where u p is the plasma frequency for non-relativistic electrons: 





Anahcn e 



(A-7) 
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For a relativistic electron plasma, where T e > m e , then [104, 105] 

47ie 2 n P m P 



2 ^ /IC - "-e _ '"-e / nonrcl\ z /a o\ 

^^^^-37;^ ) • (A " 8) 



The dimensionless function B(v) or -8(7) encodes the effects of different physical regimes. 
We are interested in various limits, specifically ion speed v ~ 1 versus ti< 1, and relativistic 
T e ^> m e versus non-relativistic T e <C ^ e plasmas. The B function generally features a 
Coulomb logarithm, which is roughly A ~ ln(g max /g m i n ), where q is the momentum transfer 
from the fast particle to a plasma electron. For example, in the non-relativistic case, g max ~ 
2m e , yv, and q min ~ hu p /v, so that A ~ ln(2m e 7t> 2 /hu p ). Then, for T e -C m e , i.e., n e -C T e 3 , 
we put 



5(7) = In 



27me ^ erff^l-^^fl + ^e-W-) 2 



(A-9) 



where u|i = 2T e /m e . The term in square brackets — > 1 for relativistic nuclides, but gives the 
appropriate reduction in losses for non-relativistic particles. 

At high energies > 1 GeV, two other loss mechanisms are important, and our treatment 
here closely follows that of [22]. Photo-pair production, e.g., py — > pe + e~ , becomes important 
for e p > 3 GeV (T/0.1 MeV) [107]. We include this using the procedure outlined in [108]. 
Energy losses due to photopion production, e.g., — > pre , become important at energies 
e p > 500 GeV (T/0.1 MeV). We include these following [109]. 

For neutrons, Coulomb scattering is of course absent, but magnetic-moment interactions 
with the ambient radiation background do act to slow the non-thermal neutrons [11,110]. 
We adopt the loss rate of [111]: 

K = = 2ne ' 9 ^ ,V, (A-10) 

at m e \m n J 

where m n is the neutron mass and g n = —1.91 is the magnetic moment of the neutron in 
units of the nuclear magneton. This expression is similar to the result of [27], except that 
the latter has a numerical prefactor 3/2 compared to that in eq. flA-lOj) and omits the factor 
(1 + v 2 ). These differences are not dramatic, and for mildly relativistic speeds these factors 
compensate to give very similar loss rates. 

A. 3 Solution of the Propagation Equation 

All the non-thermal species h evolve according to cascade equations of the form (|T6|) . and 
together these constitute a coupled set of partial differential equations (PDEs) in energy and 
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time. Because the source term includes the elastic term with Nh inside the integral, these 
equations are of integro-differential form. 

Fortunately, major simplifications occur due to the hierarchy of timescales in the problem. 
The equation is self-regulating, and the timescale to reach quasi-steady equilibrium is much 
faster than the decay time or the expansion timescale. Hence one can treat the spectrum 
as quasi-steady, so that the coupled PDEs in energy and time become coupled ordinary 
differential equations (ODEs) at each time step, a great simplification. Moreover, redshift 
and expansion effects can safely be ignored. 

Despite these considerable improvements, the equation is still of integro-differential form, 
and the integration is therefore not immediate. We follow the usual procedure to solve such 
integral equations, namely iteration. We begin with a zeroth-order approximate solution, 
in which we neglect the elastic and inelastic source terms. In this approximation, eq. (TT6l) 
becomes 

d t N h (e) = J x , h (e) - T h (e)N h (e) - d e (b h (e)N h (e)). (A-ll) 

For secondary nuclear species not present in the original decays, (i.e., all but nucleons), the 
solution of eq. ( 1A-11I) can be written in terms of the following quadrature 

1 r°° 

N h (e,t) = r— de' JxA^t) e-^>*), (A-12) 
6(e) Je 

where the exponential "optical depth" factor 

T(£ ' t)= /" rf£ "w (A " 13) 

is a measure of the average number of inelastic interactions over the time taken to lose energy 
electromagnetically from e' to e. The secondary and tertiary cascades can then be treated 
iteratively, correcting the initial omission of these terms: 

1 i>oo 

4°\e,t) = —— de' J Xth V,t)e-*W\ (A-14) 
t>h{e,t) J e 

Ng-\e,t) = — — / de' Jf(e',t)e-^), (A-15) 

': (A-16) 



OG 



iVf(e,t) = de' Jf 1) (e',t)e-^), (A-17) 

where J® is the source term from eq. ( 1A-3I) . replacing N h with N®. 

Once these distributions converge, one can insert them into equation ( fl3l) and solve for 
the hadro-dissociation rate. This iterative procedure generalizes the procedure adopted in 
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some previous studies [22]. Rather than including the exponent in the integral, they treat 
the exponential R term as a 5 function, evaluated at e*, defined as the energy at which the 
optical depth is unity: r(e*,e, t) = 1. 

We plot the resulting non-thermal neutron and proton spectra in Fig. [32] The thick solid 
black curve gives (in arbitrary units) the shape of the source spectra eQ p (e) and eQ n (e). The 
colored curves give the corresponding propagated spectra for a wide range of cosmic epochs, 
labeled by the temperature. The highest temperature, T = 2.5 MeV, corresponds to the 
lowest (red) curve, while the lowest temperature, T = 6 x 1CT 5 MeV, corresponds to the 
highest (bluegreen) curve in each plot. Note that we plot eNh(e)/Tx- (1) we use the product 
eNh(e) so that one can best assess the particle number per logarithmic interval; and (2) as 
noted above in Section [331 the propagated spectra are proportional to the decay rate T x , 
and so by dividing by this factor we remove this dependence. 

Most striking in Fig. [32] is the enormous dynamic range in the spectra amplitudes as the 
cosmic environment evolves. Indeed, the highest and lowest curves in the plots are separated 
by more than 20 orders of magnitude! This can be understood as a combination of two 
factors. First, the propagated spectra scale inversely with the cosmic baryon density, either 
directly via N^thm oc 1/n^ (eq. [20]), or indirectly via the density dependence of energy losses, 
-W/i,thick ^ ^-/b oc l/Y e ±n B (eq. [21]). This inverse density dependence leads to the steady rise 
in the spectral amplitudes at late times. Since the temperature range shown in Fig. [32] varies 
over a factor Thi/Ti ow = 4 x 10 4 , the baryon density varies by a factor (Xhi/Ti ow ) 3 = 6 x 10 13 . 
For neutrons at e < 100 MeV, the spectra are dominated by thin-target interactions and 
this factor gives the range over which these spectral amplitudes vary. 

It is important to note that while the density changes lead to enhanced spectra amplitudes 
at low temperatures and densities, the non-thermal interaction rates scale as the product of 
the non-thermal flux times the density of background targets. This compensating factor of 
density means that the overall interaction rates do not directly vary with density (though 
the spectral shapes vary with cosmic epoch). 

For neutrons at e > 100 MeV, and for protons at all energies, there is an additional 
evolutionary effect on the spectral amplitudes. Namely, these particles are dominated by 
thick-target propagation, where the spectral amplitudes are inversely proportional to energy 
losses which are themeselves dependent on the e density and thus the total lepton per baryon 
ratio Y e ±: N hjthick oc 1/6 oc l/Y e ±n B (eq. [21]). Before pair annihilation, pairs are comparable 
in number to photons, Y e ± ~ Y 1 ~ l/i] ~ 10 9 ; after annihilation is complete, the remaining 
electrons provide charge balance and thus Y e ± = Y e - ~ 1. Thus as pairs annihilate and 
energy losses are greatly reduced, the propagated spectra rise dramatically. The pair number 
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changes most when temperature drops from T ~ 0.1 MeV down to ~ 0.02 MeV, which leads 
to the large increases in the neutron spectral amplitudes at e > 100 MeV and in the proton 
amplitudes at all energies. This additional annihilation factor contributes on top of the 
baryon density effect and leads to the full variation seen in the spectral amplitudes. Note 
that since the thin-target neutrons are responsible for most of the abundance perturbations, 
and since we are interested in times > 10 3 sec, this annihilation effect is not very important 
for our results, though it is included in all of our calculations. 

Beyond the overall spectral amplitudes, the spectral shapes versus energy are set by the 
combination of the source shape and the energy dependence of the dominant loss curve. For 
example, in the case of neutrons at e < 100 MeV, the losses are dominated by inelastic 
scattering which is nearly energy-independent, and thus the nearly flat source shape eQ n (e) 
leads to a nearly flat propagated shape. On the other hand, for protons at late times and 
at e < 1 GeV, losses are dominated by Coulomb interactions, with N p oc Q p (> e)/b coa \. We 
have & C oui oc l/v oc 1/e 1 / 2 , and for this decay spectrum Q p (> e) ~ e<5 P (e) ~ constant. Thus 
we have eN p (e) oc e 3 / 2 , which is exactly the scaling behavior in Fig. [32] for < 1 GeV protons 
at late times and low temperatures. 

A. 4 Rates for Non-Thermal Reactions on Background Light Nu- 
clei 

As non-thermal particles propagate in the cosmic plasma they interact with background 
nuclides and change light element abundances. As shown in eq. the rates for these 

processes depend on the non-thermal spectra we calculate, but also on the cross sections 
for the processes. The largest such rates are between the non-thermal nucleons and the 
most abundant background species-nucleons and 4 He nuclei; these rates are important in 
determining the non-thermal spectra themselves, and are included in the propagation calcu- 
lations. Reactions with less abundant background nuclei are neglected in the non-thermal 
propagation calculations. All of the non-thermal production and rates on background nu- 
clei are evaluated as a function of cosmic time and environment, and are incorporated as 
additional channels in the usual (thermal) nucleosynthesis code. 

Table [2] lists the nuclear reactions we use in our analysis, as well as the references for the 
cross sections we adopt. We show the reactions used in the propagation calculations as well 
as in the non-thermal rate evaluations, and reactions only used for non-thermal rates. 

Several reactions have been updated since previous studies or were not present. We 
adopt the recent revised fit to 3 He(a,7) 7 Be from Cyburt and Davids [99]. The cross section 
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Figure 32: Spectra of non-thermal hadrons after propagation in the cosmic plasma; we 
plot eNh(e), which indicates the particle number per logarithmic energy interval. Results are 
shown for a large range of cosmic epochs (i.e., temperatures, densities, background abun- 
dances). Propagated spectra scale with the decay rate Tx, so we plot Nh/Tx- The decay 
injection spectra are for (mi/ 2 , m 3 / 2 , tan/3) = (120GeV, 500GeV, 50). 
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is somewhat higher than in previous assessments, leading to larger non-thermal production 
of 7 Be as well as a boosted thermal production and thus a higher standard (unperturbed) 
BBN value [6]. All else being equal, this change increase worsens the 7 Li problem, and 
also somewhat reduces the 6 Li/ 7 Li ratio at where the 7 Li level remains dominated by the 
primordial production. 

We have also added non-thermal reactions not present in our previous treatments, and as 
far as we are aware, also not present in those of other groups. Namely, we have included new 
channels for non-thermal A = 7 destruction: 7 Be(n,p) 7 Li, 7 Be(n, a) 4 He, and 7 Li(p, a) 4 He. 
Of these, by far the largest cross section is for 7 Be(n, p) 7 Li, which as a thermal rate is 
important in standard BBN. This reaction leads to A = 7 reduction when followed by 
7 Li(p, a) 4 He which has a reduced Coulomb barrier relative to 7 Be + p. As noted above 
(Section [5]) and in refs. [21,22], A = 7 destruction via 7 Be(n,p) 7 Li with thermalized decay 
neutrons is responsible for the "valley" of concordant 7 Li/H at lifetimes 10 2 — 10 3 sec as seen 
in Fig. HI It thus is reasonable to include the same process due to non-thermal neutrons. In 
fact, we find that these rates do not appreciably change the final 7 Li/H abundance. 



B Gravitino Decay Amplitudes 

In calculating the relevant gravitino decay rates, we use the interaction Lagrangian given 
in [112]: 



V2M P - 
i — 



sMp WnV^ G) F£\ (A-18) 

where ip^ is gravitino field, and x are chiral matter fields, F pa is the field strength of 
a gauge boson field, and A is the corresponding gaugino field. We denote xl = PlX an d 
Xl — xPr-, where Pi and Pr are the usual chiral projectors. The indices (C) and (G) 
denote chiral and gauge multiplets, respectively, and M P is the reduced Planck mass defined 



as Mp = l/v8vrG/v, where Gn is the Newton gravitational constant. The Feynman rules 
derived from this Lagrangian are given in Appendix B of [112]. 
We use the relations 

p L w = p l (v^xx + v 2 \x2) , p l h = p l (v* 2 xi + v; 2 x2) , 

PrW = P R (Unxi + U 21 X2) , PrH = P R (U 12 Xi + U 22 x 2 ) , (A-19) 

to transform the interaction-eigenstate charged wino W + and charged Higgsino H + to the 
mass-eigenstate charginos x+, where U and V are the unitary matrices used to diagonalize 
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the chargino mass matrix. We note that all the fields in the above relations are positively 
charged. To get the transformation relations for negatively charged conjugate states of the 
positively charged fields, one needs to exchange U and V in the above relations, e.g., 



P L W C = P L (U* n xl + E/JxxS) • (A-20) 

For neutralinos, the transformation from the interaction eigenstates Hi j2 , W 3 and B to the 
neutralino mass eigenstates Xi is given by the relations 

P L Hi = P L Nl i+2 x] , PlB = PlN^x* , PlW, = P L N* 2 x] , 

P R Hi = PrN^x] , PrB = P R N ]lX ° , PrW 3 = PrN^ , (A-21) 

where the repeated index j is summed over from 1 to 4, and the unitary matrix N diagonalizes 
the neutralino mass matrix. To obtain the mass eigenstates of the Higgs bosons, we use the 
relations [113] 

Hi = H + cos (3, Hi = H~ sin (3, 

Hi = vi + —= [Hi cos a - H% sin a + iH% sin (3) , 
V2 

Hi = v 2 + -^= (H^ sin a + H% cos a + iH® cos 0) , ( A-22) 

V2 

where the H ± are the charged Higgs fields, H^ and H® are the CP-even neutral Higgs 
fields, and H® is the CP-odd neutral Higgs field. We denote by v\ and v 2 the vacuum 
expectation values of the two neutral Higgs fields, tan/5 is defined as tan/3 = v 2 /vi, and a 
is the mixing angle which leads to the two CP-even mass eigenstates. Together with the 
relation = \ g 2 {y\ + v%), one gets vi = cos (3 and v 2 = sin (3, where g is the 

gauge coupling for SU (2) L and mw is the mass of the W boson. For squarks and sleptons, we 
do not consider mixing between generations, and we take into account only left-right mixing 
for the third generation. This simplification is sufficient within the CMSSM framework we 
consider. The relation between the mass eigenstates t^ 2 and the interaction eigenstates ti,R 
is 




(A-23) 



where the transformation matrix is unitary. There are similar interaction/mass-eigenstate 
transformations for b and f. We treat the neutrinos as in the Standard Model, i.e., as 
massless, purely left-handed neutrinos (and right-handed anti-neutrinos) only. For details of 
the above transformations, we direct the reader to [113]. 
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We now list the amplitudes for all the two-body gravitino decay channels. We denote by 



m b and m/ the masses of the final-state bosons and fermions, respectively. For G — > Xil 
and G — > g g, we have: 



\M[ 



Qml /2 Ml 



\B\ 2 (3m 2 3/2 + m 2 )( 



m 3/2 - m 



(A-24) 



where B = N' iX for the Xi final states, with N' a = N a cos 6w + N i2 sin Q Wl and 6w is the weak 
mixing angle, and B = 1 for the g final state. For G ^ xt W T and G — > we have: 



12m 3/2 



j— T 1 12 (CD + C*D*) ml /2 mlm f 



2^ 



+ (1^1 +1^1) 3 m 3/2 - m 3/2 ( 5 m / + m t) + m 3/2 ( m / - m fc) + ( m / - m 6) 

c 2 



24 m\ j2 mlMl 



2^2 



m 3/2 - 2 m 3/2 ( m / - 5 m 2 ) + (mj - m 2 ) 



X 



+ 



X 



2 (E F + E*F*) m 3/2 m f + (\E\ 2 + \F\ 2 ) (m 2 3/2 + m) - m 2 b ) 



G 



6m 3/2 Mp { 



(DE* + ED* + CF* + FC*) 

q4 i 2/2, 2\ , / 2 2\ 2 

-2m 3/2 + m 3/2 [m f + m b ) + [m f - m b ) 
3(CE + DF + C*E* + D*F*) m 3/2 m f (m 2 3/2 - m) + m 2 ) j, 



(A-25) 



where C = Vji, D = Uji, E = —viUj 2 , F = —v 2 Vj 2 and G = g for the xf nna l states; 
and C = D = A^ 2 , with N[ 2 = -A a sin0 w + N i2 cos6 w , E = F = (-v^ + v 2 N i4 ) 
and G = gj cos#vk for the Xi final states. In the above expression, the first and second 
lines come from the gauge part of eq. (IA-180 (the second line), the third and fourth lines 
come from the matter part of eq. (1A-18j) (the first line), with the gauge bosons coming 
from the covariant derivative and the Higgs field <fi taking its vacuum expectation value (see 
eq. ( 1A-22I) ). and the fifth and sixth lines come from the interference of these two parts. For 
G - / /, G - and G - xl K 



: " a ° 2 3 , we have: 



Ylm\ l2 Ml 



(\H\ 2 + \K\ 2 ) (m 2 3/2 + m 2 f -m 2 b )-2(HK + H*K*) m 3/2 m f 



m i/2 - 2 m 2 3/2 (mj + m 2 b ) + (m) - mf)' 



(A-26) 



where H = K = (A i4 sina + A i3 cosa) for the H® final state; H = K = ^(A^cosa— 
N i3 sin a) for the H 2 final state; H = —K = 4| (Na cos (3 + N i3 sin /3) for the H 3 final state; 
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H = Uj2 sin (3 and K = Vj2 cos f3 for the H final state; H = U~. R and K = Uj kL for the 
tk final state (and similarly for bk and fk with their corresponding mixing matrices and 
Uf, respectively); H = and = 1 for all the fi and z> final states; and H = 1 and = 
for all the final states. To get the partial widths, one need to multiply (1A-26j) by the 
phase-space factor: 



16 * m l/2 



711 



3/2 



-(m/-m b ) ml /2 - (rrif + m b y 



(A-27) 



where N c is the color factor (3 for qq channels, 8 for the gg channel, and 1 otherwise). 
The amplitude for the three-body decay G 

g 2 sin 2 OwQl 



Xi QQ is 



\M\< 



3 s 2 m 2 /2 M 2 



Asm n ml /2 (N>? + N>1 2 ) (s + 2m 2 q ) 



+ \K\ 2 



3m 3/2 | s 



2ml) +^3/ 2 (-3s(s + 2t) +m|u (s-10m 2 )) 

+ m l/2 (6 s m\ + 2 m 2 ( - 3 s (s + 2 1) + 2 s m|o + m\o) 

+ s(s 2 + 8st + 6t 2 -4m|o (s + 2 1) + 3 m|o)) - (s - m|.) (2 s mj - 2 



x 



(2 + s m 2 - ml) + s (s 2 + 2 st + 2 1 2 - 2 m 2 (s + t) + mL) 



(A-28) 



where s and t are the invariant masses of the qq and systems, respectively, m^o is the 
mass of the neutralino where the LSP corresponds to i — 1, and Q g = 2/3 or —1/3 for 
the corresponding quarks. The differential decay rate is 

dT N r 



ds dt 



256 7T 3 m^ 2 



\M\' 



(A-29) 



These three-body and two-body analytic expressions were also given in [24] and [42]. 

For G — > x < iW + W~, there are contributions from four generically different diagrams: 
contact, j/Z exchange, xf exchange and Hf/H 2 exchange diagrams. The amplitude is 
too long to be listed here, so we just write the matrix elements for each of these diagrams 
(suppressing the polarization indices): 



iM 



contact 



4M, 



■gu(q') (P R N l2 + P L N* 2 ) -f [l P , la] Mp) e°*(k') 



iM 



G-*x i -y*-+Xi w+w- 



AM P (k + k 



— g^g sin 8 W u(q') {P R N' a + P L N'*) -f [t + H' , 1*} Mp) 



x 



(2 k + k') a gr3p + (k' - k) g pu - (k + 2 k') g af3 e»*(k) e°*(k') 
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iM G ^o z *^o w + w - 



AM P (k + k') 2 - m| + im z T z 
x u(q') | (P R N' l2 + P L N'*) Y [t + ti\ la 
V2g 



9 p — ) 9 cos 6 W 



Tflr 



X 



COS0 



w 



Pr (-V!N i3 + v 2 N l4 ) + P L (- Vl N* 3 + v 2 N* A ) Yla MP) 



(2 k + k') a g Pp + (k' - k) p g ap — (k + 2 k') p g^ e?*(k) e°*{k') , 



gu(q') la [P R 0« + P L OL] 



4M P (p - k) 2 - m 2 ± + im-±T-± 

Xj x j x j 

■(f-k , + m^) | {p R v n + PlU*,) rh P , *] 

+ 2g {P R v 2 V j2 + P LVl W j2 ) Yip] MP) t P *(k) e°*(k') , 



iM 



gu{q') lp [P L Of j * + P R 0^] 



4Mp (p - k') 2 - m 2 ± + im~±T~± 

Xj X j Xj 

■(f-k" + m-±) | (P^* + P R U n ) >f[ 1<r , k"} 

+ 2g (P L v 2 V* 2 + P RVl U j2 ) >f la \ ^(p) e p *(k) e°*{k') , 



iM 



g-+x1ho*-+x'}w+w- 



■ g m w cos ((3 - a) g pa 



2M P (k + k') 2 - m 2 H o + im H oY H o 
x u(q ) [Pr (sin aN iA + cos aN i3 ) — P L (sin aN* A + cos aN* 3 ) ] 7 M 
• ^ + k")Mp)^*(k)e' 7 *(k'), 



i -^G^x i H°*^x° i W+W- 



g m w sin {(3 - a) g pa 



2M P (k + k') 2 - m^o + im H oT H o 
x w(g') [Pr (cos aN i4 - sin ceA^) - P L (cos aiV* 4 - sin aN* 3 ) ] Y 
■ ^ + k")Mp)^*(k)^*(k'), (A-30) 
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where ip^{p) represents the decaying gravitino with momentum p, u{q') represents the pro- 
duced neutralino with momentum q', e p *(k) and e a *(k') are the polarization four- vectors for 
the two W bosons with momenta k and k', respectively, g a/3 is the flat-space Lorentz metric 
tensor, the T factors in the propagators are the total widths of the exchanged particles and 
Og = NaVji - -^N i4 V* 2 and Og = N^U fl + -^N* 3 U j2 . The Feynman rules for the vertices 
W + W~H® 2 and xfXiW T are given in [113,114]. The sum of the polarization states for 
gravitino is 

E «w = -<,+ ™ 3/2 ) (,„ - (7, + (f- M (7„ + ^) , 

and the phase space is the same as given in eq. (IA-291) . 

A remark is in order for the process G — > xl W + W~ . Because of the 1/mw factors in 
the longitudinal polarization states of the W boson, one might worry that this process might 
have bad high-energy behavior. Correspondingly, if one restored the electroweak symmetry to 
make mw vanish, then this process would diverge if there were still some terms proportional 
to negative powers of my/ in the final result. However, upon expanding the parameters in 
terms of mw, that is, expanding the elements of chargino and neutralino mixing matrices 
N, U, V, the masses m^o, m~±, and the angle a etc. (for the mw dependence of these 

A-i Xj 

parameters, see, for example [115]), in the total amplitude, we have verified that all terms 
with negative powers of mw cancel. As a result, this process has good high energy behavior. 
This is a highly non-trivial check on our calculation of the process G — > x® W + W~ , verifying 
the relative magnitudes, signs and phases of all the individual contributions to the decay 
amplitudes. 

We conclude by commenting on the differences between the amplitudes obtained and 
used this paper and those given [24] and [42], which are based on the gravitino interaction 
Lagrangian given in eq. (4.58) of [37]. This Lagrangian is different from that given in 
eq. (2.82) of [112], which we follow here. The latter reference uses the transformation from the 
two-component gravitino Lagrangian given in Wess and Bagger [116] to the four-component 
Lagrangian which, as pointed out at the footnote on page 205 of [114], requires a factor 
of i when one defines the Majorana spinors. This factor results in a sign difference when 
one writes down the Feynman rules using the matter part of the Lagrangian. We also note 
that our result for gravitino — ► Z + x differs from [42], in that it includes the Higgsino 
contribution. 



64 



References 



[1] J. Dunkley et al. [WMAP Collaboration], Astrophys. J. accepted [ arXiv:0803.0586l astro- 

ph]. 

[2] R. H. Cyburt, B. D. Fields and K. A. Olive, New Astron. 6 (2001) 215 
[arXiv:astro-ph/01021~79) . 

[3] R. H. Cyburt, B. D. Fields and K. A. Olive, Astropart. Phys. 17 (2002) 87 
|arXiv:astro-ph/0105397| . 

[4] R. H. Cyburt, B. D. Fields and K. A. Olive, Phys. Lett. B567, 227 (2003); 
A. Coc, E. Vangioni-Flam, P. Descouvemont, A. Adahchour and C. Angulo, Astro- 



[5] 
[6] 

[7 
[8 
[9 
[10 

[11 
[12 



[13 
[14 



phys. J. 600 (2004) 544 [arXiv:astro-ph/ 0309480 1 ; A. Cuoco, F. Iocco, G. Mangano, 
G. Miele, O. Pisanti and P. D. Serpico, Int. J. Mod. Phys. A 19 (2004) 4431 
|arX iv:astro-ph/0307213| ; B.D. Fields and S. Sarkar in: C. Amsler et al. [Particle Data 
Group], Phys. Lett. B 667, 1 (2008); P. Descouvemont, A. Adahchour, C. Angulo, 
A. Coc and E. Vangioni-Flam, ADNDT 88 (2004) 203 |arXiv:astro-ph/040710 1|. 

R. H. Cyburt, Phys. Rev. D 70, 023505 (2004) [arXiv:astro-ph/0401091] . 

R. H. Cyburt, B. D. Fields and K. A. Olive, JCAP 0811 (2008) 012. [arXiv:0808.2"8l"8l 
[astro-ph]]. 

D. Lindley, Astrophys. J. 294 (1985) 1. 

J. R. Ellis, D. V. Nanopoulos and S. Sarkar, Nucl. Phys. B 259 (1985) 175. 

D. Lindley, Phys. Lett. B 171 (1986) 235. 

R. J. Scherrer and M. S. Turner, Astrophys. J. 331 (1988) 19. 

M. H. Reno and D. Seckel, Phys. Rev. D 37 (1988) 3441. 

S. Dimopoulos, R. Esmailzadeh, L. J. Hall and G. D. Starkman, Astrophys. J. 330, 545 
(1988); S. Dimopoulos, R. Esmailzadeh, L. J. Hall and G. D. Starkman, Nucl. Phys. B 
311 (1989) 699. 

J. Ellis et al, Nucl. Phys. B 337, 399 (1992). 

M. Kawasaki and T. Moroi, Prog. Theor. Phys. 93 (1995) 879 |arXiv:hep-ph~79~403364 ] . 



65 



[15 
[16 

[it; 

[18 

[19 
[20 

[21 

[22 

[23 

[24 

[25 

[26 

[27; 

[28 
[29 



M. Kawasaki and T. Moroi, Astrophys. J. 452, 506 (1995). 

E. Holtmann, M. Kawasaki, K. Kohri and T. Moroi, Phys. Rev. D 60, 023506 (1999) 
[arXiv:hep-ph/9805405| . 



K. Jedamzik, Phys. Rev. Lett. 84, 3248 (2000) |arXiv:astro-ph/9909445|. 

M. Kawasaki, K. Kohri and T. Moroi, Phys. Rev. D 63 (2001) 103502 
arXiv:hep-ph/0012279| . 



K. Kohri, Phys. Rev. D 64 (2001) 043515 |arXiv:astro-ph/0103411] . 

R. H. Cyburt, J. R. Ellis, B. D. Fields and K. A. Olive, Phys. Rev. D 67 (2003) 103521 
arXiv:astro-ph/0211258~| . 



K. Jedamzik, Phys. Rev. D 70 (2004) 063524 |arXiv:astro-ph/0402344 j[ ; K. Jedamzik, 
Phys. Rev. D 70 (2004) 083510 [arXiv:astro-ph/0405583| . 

M. Kawasaki, K. Kohri and T. Moroi, Phys. Lett. B 625 (2005) 7 
arXiv:astro-ph/0402490| ; Phys. Rev. D 71 (2005) 083502 |arXiv:astro-ph/0408426| . 

J. R. Ellis, K. A. Olive and E. Vangioni, Phys. Lett. B 619, 30 (2005) 
|arXiv:astro-ph/0503023] . 

K. Kohri, T. Moroi and A. Yotsuyanagi, Phys. Rev. D 73, 123511 (2006) 



arXiv:hep-ph/0507 245| . 

D. G. Cerdeno, K. Y. Choi, K. Jedamzik, L. Roszkowski and R. Ruiz de Austri, JCAP 



0606, 005 (2006) |arXiv:hep-ph705 09275|. 

K. Jedamzik, K. Y. Choi, L. Roszkowski and R. Ruiz de Austri, JCAP 0607, 007 (2006) 
arXiv:hep-ph/0512044|. 



K. Jedamzik, Phys. Rev. D 74, 103509 (2006) |arXiv:hep-ph/0604251 



F. D. Steffen, JCAP 0609, 001 (2006) |arXiv:hep-ph7 0605306|. 

R. H. Cyburt, J. R. Ellis, B. D. Fields, K. A. Olive and V. C. Spanos, JCAP 0611, 014 
(2006) |arXiv:astro-ph/0608562|. 



[30] A. Bouquet and P. Salati, Nucl. Phys. B 284, 557 (1987); J. E. Kim, B. Kyae and 
J. D. Park, |arXiv:hep-ph/9810503| 

66 



[31] J. Ellis, J.S. Hagelin, D.V. Nanopoulos, K.A. Olive and M. Srednicki, Nucl. Phys. B238, 
453 (1984). 



[32 
[33 
[34 
[35 
[36 
[37; 
[38 

[39 

[40 
[41 



[42 
[43 
[44; 



[45 



D. V. Nanopoulos, K. A. Olive and M. Srednicki, Phys. Lett. B 127 (1983) 30. 

J. R. Ellis, J. E. Kim and D. V. Nanopoulos, Phys. Lett. B 145 (1984) 181. 

R. Juszkiewicz, J. Silk and A. Stebbins, Phys. Lett. B 158 (1985) 463. 

M. Kawasaki and K. Sato, Phys. Lett. B 189 (1987) 23. 

T. Moroi, H. Murayama and M. Yamaguchi, Phys. Lett. B 303 (1993) 289. 

T. Moroi, |arXiv:hep-ph/9503210| 

J. R. Ellis, D. V. Nanopoulos, K. A. Olive and S. J. Rey, Astropart. Phys. 4 (1996) 371 
|arXiv:hep-ph/9505438l . 

M. Bolz, A. Brandenburg and W. Buchmuller, Nucl. Phys. B 606 (2001) 518 
|arXiv:hep-ph/0012052| . 

J. Pradler and F. D. Steffen, Phys. Lett. B 648, 224 (2007) | arXiv:hep-ph/0612291] . 

J. L. Feng, A. Rajaraman and F. Takayama, Phys. Rev. D 68 (2003) 063504 
[arXiv:hep-ph/0306024 | ; J. L. Feng, S. F. Su and F. Takayama, Phys. Rev. D 70 (2004) 
063514 [arXiv:hep-ph/0404198] . 

J. L. Feng, S. Su and F. Takayama, Phys. Rev. D 70 (2004) 075019 
|arXiv:hep-ph/0404231| . 

J. R. Ellis, K. A. Olive, Y. Santoso and V. C. Spanos, Phys. Lett. B 588 (2004) 7 
|arXiv:hep-ph/0312262| . 

J. R. Ellis, K. A. Olive, Y. Santoso and V. C. Spanos, Phys. Lett. B 573 (2003) 162 
[arXiv:hep-ph/0305212 |; J. R. Ellis, K. A. Olive, Y. Santoso and V. C. Spanos, Phys. 
Rev. D 70 (2004) 055005 [arXiv:hep-ph/0405110l . 

J. R. Ellis, K. A. Olive, Y. Santoso and V. C. Spanos, Phys. Lett. B 565, 176 (2003) 
|arXiv:hep-ph/0303043l . 



67 



[46] H. Baer and C. Balazs, JCAP 0305, 006 (2003) |arX iv:hep-ph/03031l4 |; A. B. La- 
hanas and D. V. Nanopoulos, Phys. Lett. B 568, 55 (2003) [arXiv:hep-ph/0303130| ; 
U. Chattopadhyay, A. Corsetti and P. Nath, Phys. Rev. D 68, 035005 (2003) 
[arXiv:hep-ph/0303201 |; C. Munoz, Int. J. Mod. Phys. A 19, 3093 (2004) 
[arXiv:hep-ph/0309346| ; R. Arnowitt, B. Dutta and B. Hu, |arXiv:hep-ph/0310103 

[47] M. Battaglia et al, Eur. Phys. J. C 22, 535 (2001) | arXiv:hep-ph/0106204| ; 
M. Battaglia, A. De Roeck, J. R. Ellis, F. Gianotti, K. A. Olive and L. Pape, Eur. 
Phys. J. C 33, 273 (2004) |arXiv:hep-ph/0306219| . 

[48] R. I. Epstein, J. M. Lattimer and D. N. Schramm, Nature 263, 198 (1976); T. Pro- 
danovic and B. D. Fields, Astrophys. J. 597, 48 (2003) larXiv:astro-ph/0307183l . 

[49] J. L. Linsky et al, Astrophys. J. 647, 1106 (2006) |arXiv:astro-ph/0608308] . 

[50] E. Vangioni-Flam, K. A. Olive and N. Prantzos, Astrophys. J. 427, 618 (1994) 
[arXiv: astro-ph /93 1 002 1] ; S. Scully, M. Casse, K. A. Olive and E. Vangioni-Flam, As- 
trophys. J. 476, 521 (1997) |arXiv:astro-ph79 607106|; G. Stei gman, D. Romano and 
M. Tosi, Mon. Not. Roy. Astron. Soc. 378, 576 (2007) [arXiv:astro-ph/0703682| ; T. Pro- 
danovic and B. D. Fields, JCAP 0809, 003 (2008) |arXiv:0804.3095l [astro-ph]]. 

[51] S. Buries and D. Tytler, Astrophys. J. 499, 699 (1998) |arXiv:astro-ph/9712108] . 

[52] S. Buries and D. Tytler, Astrophys. J. 507, 732 (1998) |arXiv:astro-ph/9712109] . 

[53] J. M. O'Meara, D. Tytler, D. Kirkman, N. Suzuki, J. X. Prochaska, D. Lubin and 
A. M. Wolfe, Astrophys. J. 552, 718 (2001) [arXiv: astro-ph/00 1 1 1 79| . 

[54] M. Pettini and D. V. Bowen, Astrophys. J. 560, 41 (2001) |a rXiv:astro-ph/0104474| . 

[55] D. Kirkman, D. Tytler, N. Suzuki, J. M. O'Meara and D. Lubin, Astrophys. J. Suppl. 
149, 1 (2003) [arXiv:astro-ph/0302006] . 

[56] J. M. O'Meara, S. Buries, J. X. Prochaska, G. E. Prochter, R. A. Bernstein and 
K. M. Burgess, Astrophys. J. 649, L61 (2006) |arXiv: astro-ph/0608302| . 

[57] M. Pettini, B. J. Zych, M. T. Murphy, A. Lewis and C. C. Steidel, MNRAS 391, 1499 
(2008) |arXiv:0805.0594l [astro-ph]]. 

[58] Bania, T. M., Rood, R. T., & Balser, D. S. 2002, Nature, 415, 54. 



68 



[59 

[60 

[61 

[62 

[63 
[64; 
[65 
[66 

[67; 

[68 
[69 

[70 
[71 
[72 

[73 
[74 



E. Vangioni-Flam, K. A. Olive, B. D. Fields and M. Casse, Astrophys. J. 585, 611 
(2003) |arXiv:astro-ph/0207583| . 

B. D. Fields, K. A. Olive, J. Silk, M. Casse and E. Vangioni-Flam, Astrophys. J. 563, 



653 (2001) |arXiv:astro-ph7 0107389|. 

G. Sigl, K. Jedamzik, D. N. Schramm and V. S. Berezinsky, Phys. Rev. D 52 (1995) 
6682 |arXiv:astro-ph/9503094| . 

M. Peimbert, A. Peimbert and M.T. Ruiz, Astrophys. J. 541, 688 (2000); A. Peimbert, 
M. Peimbert and V. Luridiana, Astrophys. J. 565, 668 (2002). 

Y. I. Izotov and T. X. Thuan, Astrophys. J. 500, 188 (1998). 

K.A. Olive, and E. Skillman, New Astronomy, 6, 119 (2001). 

Y. I. Izotov and T. X. Thuan, Astrophys. J. 602, 200 (2004) [arXiv: astro-ph/031042 1| . 



K. A. Olive and E. D. Skillman, Astrophys. J. 617, 29 (2004) |arXiv:astro-ph/0405588 



M. Fukugita and M. Kawasaki, Astrophys. J. 646, 691 (2006) |arXiv:astro^p h/0603334|. 
F. Spite, M. Spite, Astronomy & Astrophysics, 115 (1992) 357. 

S. G. Ryan, T. C. Beers, K. A. Olive, B. D. Fields, and J. E. Norris Astrophys. J. Lett. 
530 (2000) L57 [arXiv: astro-ph/990521 1] . 

P. Bonifacio et al, Astron. Astrophys., 390, 91 (2002). |arXiv:astro-ph/0204332] . 
L. Pasquini and P. Molaro, Astron. Astrophys. 307, 761 (1996). 

F. Thevenin, C. Charbonnel, J. A. d. Pacheco, T. P. Idiart, G. Jasniewicz, P. de Laverny 
and B. Plez, Astron. Astrophys. 373, 905 (2001) | arXiv:astro-ph/0105166| . 

P. Bonifacio, Astron. Astrophys. 395, 515 (2002) |arXiv:astro-ph/0209434l . 

A. Hosford, S. G. Ryan, A. E. G. Perez, J. E. Norris and K. A. Olive, Astron. Astrophys. 
493, 601 (2009) |arXiv:0811.2506l [astro-ph]]. 



[75] S. Vauclair and C. Charbonnel, Astrophys. J. 502, 372 (1998); M. H. Pinson- 
neault, T. P. Walker, G. Steigman and V. K. Narayanan, Astrophys. J. 527, 180 
(2002) [arXiv:astro-ph/9803073|; M. H. Pinsonneault, G. Steigman, T. P. Walker and 



69 



V. K. Narayanan, Astrophys. J. 574, 398 (2002) |arXiv:astro-ph/0105439| ; O. Richard, 
G. Michaud and J. Richer, Astrophys. J. 619, 538 (2005) [arXiv: astro-ph/0409672| ; 
A. J. Korn et al., Nature 442 (2006) 657 | arXiv:astro-ph/0608201] . 

[76] V.V. Smith, D.L. Lambert, and P.E. Nissen, Astrophys. J. 408, 262 (1993); Astrophys. 
J. 506, 405 (1998); L.M. Hobbs and J.A. Thorburn, Astrophys. J. Lett., 428, L25 
(1994); Astrophys. J. 491, 772 (1997); R. Cayrel, M. Spite, F. Spite, E. Vangioni-Flam, 
M. Casse, and J. Audouze, Astron. Astrophys. 343, 923 (1999). 



77 



79 

80 
1 

82 
83 



84 



85 



86 



M. Asplund, D. L. Lambert, P. E. Nissen, F. Primas and V. V. Smith, Astrophys. J. 
644, 229 (2006) |arXiv:astro-ph/0510636| . 

D. Thomas, D. N. Schramm, K. A. Olive and B. D. Fields, Astrophys. J. 406, 569 



(1993) [arXiv:astro-p h/9206002|. 

E. Vangioni-Flam, M. Casse, R. Cayrel, J. Audouze, M. Spite, and F. Spite, New 
Astronomy, 4, 245 (1999) [ arXiv:astro-ph/9811327] . 

R. Cayrel et al, Astronomy & Astrophys. 473, L37 (2007). 

G. Steigman, B. D. Fields, K. A. Olive, D. N. Schramm and T. P. Walker, Astrophys. 
J. 415, L35 (1993). 

B.D.Fields and K.A. Olive, New Astronomy, 4, 255 (1999) [arXiv: astro-ph/98 1 1 183| . 

M. Kusakabe, T. Kajino and G. J. Mathews, Phys. Rev. D 74, 023526 (2006) 
|arXiv:astro-ph/0605255|; M. Kusakabe, T. Kajino, R. N. Boyd, T. Yoshida and 
G. J. Mathews, Phys. Rev. D 76, 121302 (2007) [a rXiv:0711.3854l [astro-ph]]: M. Kusak- 
abe, T. Kajino, R. N. Boyd, T. Yoshida and G. J. Mathews. larXiv:0711.3 858 [astro-ph]. 

K. Jedamzik, Phys. Rev. D 77, 063524 (2008) (arXiv: 0707.20701 [astro-ph] ] : K. Jedamzik, 
JCAP 0803, 008 (2008) |arXiv:0710.5153 [hep-ph]]. 

D. Cumberbatch, K. Ichikawa, M. Kawasaki, K. Kohri, J. Silk and G. D. Starkman, 
Phys. Rev. D 76, 123005 (2007) [arXiv: 0708 .00951 [astro-ph]]. 

S. Bailly, K. Jedamzik and G. Moultaka. la?Xiv:0812.0788 [hep-ph]; S. Bailly, K. Y. Choi, 
K. Jedamzik and L. Roszkowski, JHEP 0905, 103 (2009) |arXiv:0903.3"974l [hep-ph]]. 



70 



[87] E. Rollinde, E. Vangioni-Flam and K. A. Olive, Astrophys. J. 627, 666 (2005) 
|arXiv:astro-ph/0412426| ; E. Rollinde, E. Van gioni and K. A. Olive, Astrophys. J. 651, 
658 (2006) |arXiv:astro-ph/0605633] ; E. Rollinde, D. Maurin, E. Vangioni, K. A. Olive 
and S. Inoue, Astrophys. J. 673, 676 (2008) [ arXiv:0707.2"086l [astro-ph]]; T. Pro- 
danovic and B. D. Fields, Phys. Rev. D 76, 083003 (2007) |arXiv:0709.3"300l [astro-ph]]; 
M. Kusakabe. ErXiv:0803.340l [astro-ph]. 

[88] M. Pospelov, Phys. Rev. Lett. 98, 231301 (2007) | arXiv:hep-ph/0605215| ; 
C. Bird, K. Koopmans and M. Pospelov, Phys. Rev. D 78, 083010 (2008) 
|arXiv:hep-ph/07030'96| . 

[89] K. Kohri and F. Takayama, Phys. Rev. D 76, 063507 (2007) |arXiv:hep-ph/0605243| . 

[90] M. Kaplinghat and A. Rajaraman, Phys. Rev. D 74, 103004 (2006) 
|arXiv:astro-ph/0606209| . 

[91] K. Hamaguchi, T. Hatsuda, M. Kamimura, Y. Kino and T. T. Yanagida, Phys. Lett. B 
650, 268 (2007) larXiv:hep-ph/0702274| . 

[92] T. Jittoh, K. Kohri, M. Koike, J. Sato, T. Shimomura and M. Yamanaka, Phys. Rev. 
D 76, 125023 (2007) [ arXiv:0704.2914l [hep-ph]]. 

[93] J. Pradler and F. D. Steffen, Phys. Lett. B 666, 181 (2008) jarXiv:0710.2213l [hep-ph]]: 
J. Pradler and F. D. Steffen, Eur. Phys. J. C 56, 287 (2008) [arXiv:0710.4548l [hep-ph]]; 
F. D. Steffen, Phys. Lett. B 669, 74 (2008) [arXiv: 0806 .32661 [hep-ph]]. 

[94] M. Pospelov. la7Xiv:0712.0647l [hep-ph]. 

[95] M. Pospelov, J. Pradler and F. D. Steffen, JCAP 0811, 020 (2008) [ arXiv:0807.4"287l 
[hep-ph]]. 

[96] J. L. Diaz-Cruz, J. R. Ellis, K. A. Olive and Y. Santoso, JHEP 0705, 003 (2007) 
[arXiv:hep-ph/0701229| ; K. Kohri and Y. Santoso, arXiv:0811 .11191 [hep-ph]. 

[97] M. Meneguzzi, J. Audouze, and H. Reeves, Astronomy & Astrophys., 15 337 (1971); 
R. Cowsik and L. W. Wilson, In Denver 1973, Cosmic Ray Conference Vol.1, Denver 
1973, 500-505; R. Ramaty and R. E. Lingenfelter, in Solar gamma-, X-, and EUV 
radiation, ed. S. R. Kane, 363 (1975). 



71 



[98] T. Sjostrand, P. Eden, C. Friberg, L. Lonnblad, G. Miu, S. Mrenna and E. Norrbin, 
Comput. Phys. Commun. 135 (2001) 238 | arXiv:hep-ph/0010017| . 

[99] R. H. Cyburt and B. Davids, Phys. Rev. C 78, 064614 (2008) [arXiv: 0809 .32401 [nucl-ex]] . 



[100 
[101 

[102 
[103 
[104 
[105 
[106 
[107 
[108 
[109 
[110 
[111 
[112 
[113 
[114 
[115 
[116 



J. P. Meyer, Astron. & Astrophys. Suppl. 7, 417 (1972). 

S. Ando, R. H. Cyburt, S. W. Hong and C. H. Hyun, Phys. Rev. C 74, 025809 (2006) 



|arXiv:nucl-th705Tl 074| . 



P. Mohr et al, Phys. Rev. C 50, 1543 (1994). 

W. A. Fowler and F. Hoyle, Astrophys. J. Suppl. 9, 201 (1964). 

V. N. Tystovich, Sov. Phys. JETP 15, 561 (1962) [Zh. Eksp. Teor. Fiz. 42, 803 (1962)]. 
R. J. Gould, Physica 58 (1972) 379. 

M. Inokuti, Y. Itikawa and J. E. Turner, Rev. Mod. Phys. 50, 23 (1978). 

G. R. Blumenthal, Phys. Rev. D 1, 1596 (1970). 

M. J. Chodorowski, A. A. Zdziarski, and M. Sikora, Astrophys. J. 400, 181 (1992) 

V. S. Berezinsky and S. I. Grigor'eva, Astron. Astrophys. 199, 1 (1988). 

R. J. Gould, Astrophys. J. 417, 12 (1993). 

J. E. Turner et al, Phys. Rev. B 8 4057 (1973). 

J. Pradler. la7Xiv:0708.2786 j [hep-ph]. 

J. F. Gunion and H. E. Haber, Nucl. Phys. B 272 (1986) 1. 

H. E. Haber and G. L. Kane, Phys. Rept. 117 (1985) 75. 

J. F. Gunion and H. E. Haber, Phys. Rev. D 37 (1988) 2515. 

J. Wess and J. Bagger, Supersymmetry and supergravity, (Princeton University Press, 
Princeton, USA, 1992). 



72 



